#===========================================================================================================
# Date: March 18, 2026
# European Journal of Political Economy EJPE
# Authors: Ingrid Mauerer, Andrea Pili, Socorro Puy
# Title: What Drives Pro-EU Preferences? Territorial Identities, Immigration and Redistribution Attitudes
#===========================================================================================================

#==================
# Packages
#==================
library(ggplot2)
library(VGAM)
library(lme4)

#==================
# Load data
#==================
list2env(readRDS("Data.RDS"), envir = .GlobalEnv)


#--------------------------------------------------------
# Figure 2: Support for EU Membership across Countries
#--------------------------------------------------------

pdf(file="Figure_2.pdf", width=11, height=5)

# Matrix
country_eus <- c(rep("Germany",2),rep("France",2),rep("Italy",2),rep("Spain",2),rep("Poland",2),rep("Romania",2),rep("Netherlands",2),
                 rep("Sweden",2),rep("Hungary",2),rep("Austria",2),rep("Bulgaria",2),rep("Slovakia",2),rep("Latvia",2),rep("Estonia",2))

eus  <- rep(c("Good", "Bad"),14)

eus_value <- c(
  # DE  
  prop.table(table(DE$eu_support))[2],prop.table(table(DE$eu_support))[1],
  # FR
  prop.table(table(FR$eu_support))[2],prop.table(table(FR$eu_support))[1],
  # IT
  prop.table(table(IT$eu_support))[2],prop.table(table(IT$eu_support))[1],
  # ES 
  prop.table(table(ES$eu_support))[2],prop.table(table(ES$eu_support))[1],
  # PL 
  prop.table(table(PL$eu_support))[2],prop.table(table(PL$eu_support))[1],
  # RO 
  prop.table(table(RO$eu_support))[2],prop.table(table(RO$eu_support))[1],
  # NL 
  prop.table(table(NL$eu_support))[2],prop.table(table(NL$eu_support))[1],
  # SE 
  prop.table(table(SE$eu_support))[2],prop.table(table(SE$eu_support))[1],
  # HU 
  prop.table(table(HU$eu_support))[2],prop.table(table(HU$eu_support))[1],
  # AT 
  prop.table(table(AT$eu_support))[2],prop.table(table(AT$eu_support))[1],
  # BG 
  prop.table(table(BG$eu_support))[2],prop.table(table(BG$eu_support))[1],
  # SK
  prop.table(table(SK$eu_support))[2],prop.table(table(SK$eu_support))[1],
  # LV 
  prop.table(table(LV$eu_support))[2],prop.table(table(LV$eu_support))[1],
  # EE 
  prop.table(table(EE$eu_support))[2],prop.table(table(EE$eu_support))[1])


# save as data frame and adjust order
eu_support_data <-data.frame(country_eus, eus, eus_value)
eu_support_data$country_eus <- factor(eu_support_data$country_eus, c("Germany","France","Italy","Spain","Poland","Romania","Netherlands","Sweden","Hungary","Austria","Bulgaria","Slovakia","Latvia","Estonia"))
eu_support_data$eus <- factor(eu_support_data$eus, c("Bad","Good"))
eu_support_data


# export
ggplot(eu_support_data, aes(fill=eus, y=eus_value, x=country_eus)) + 
  geom_bar(position="stack", stat="identity", color="black") +
  theme_bw()+theme(legend.position = "top")  +
  labs(y = "Percent", x="",fill="") +
  theme(axis.text=element_text(size = 12, angle = 45, hjust = 1), axis.title=element_text(size=12),
        legend.text=element_text(size=12), legend.title=element_text(size=14)) +
  scale_fill_discrete(breaks=c("Good", "Bad"),
                      type=c("gray","black")
  )


dev.off()


#------------------------------------------------------------------------------------
# Table 2: Descriptive Statistics on Ideological and Territorial Identity Variables
#------------------------------------------------------------------------------------

Table2 <- matrix(NA, 14, ncol = 14)
colnames(Table2) <- c("Immigration","sd","Redistribution","sd", "Cor(x,y)", 
                      "National Identity","sd","Regional Identity","sd", "European Identity","sd",
                      "Cor(n,r)", "Cor(n,e)", "Cor(r,e)")
rownames(Table2) <- c("Germany", "France", "Italy", "Spain", "Poland", "Romania", "Netherlands", "Sweden", "Hungary", "Austria", "Bulgaria", "Slovakia", "Latvia", "Estonia")

# DE
Table2[1,] <-round(c(mean(DE$pro_imm),sd(DE$pro_imm), mean(DE$pro_red),sd(DE$pro_red), cor(DE$pro_imm,DE$pro_red),
                     mean(DE$nat_id),sd(DE$nat_id), mean(DE$reg_id),sd(DE$reg_id),mean(DE$eur_id), sd(DE$eur_id),
                     cor(DE$nat_id,DE$reg_id), cor(DE$nat_id,DE$eur_id), cor(DE$reg_id,DE$eur_id)), digits=2)
# FR
Table2[2,]  <-round(c(mean(FR$pro_imm),sd(FR$pro_imm),mean(FR$pro_red),sd(FR$pro_red), cor(FR$pro_imm,FR$pro_red),
                      mean(FR$nat_id),sd(FR$nat_id), mean(FR$reg_id),sd(FR$reg_id), mean(FR$eur_id),sd(FR$eur_id),
                      cor(FR$nat_id,FR$reg_id), cor(FR$nat_id,FR$eur_id), cor(FR$reg_id,FR$eur_id)), digits=2)
# IT
Table2[3,]  <-round(c(mean(IT$pro_imm),sd(IT$pro_imm),  mean(IT$pro_red),sd(IT$pro_red), cor(IT$pro_imm,IT$pro_red),
                      mean(IT$nat_id),sd(IT$nat_id), mean(IT$reg_id),sd(IT$reg_id), mean(IT$eur_id),sd(IT$eur_id),
                      cor(IT$nat_id,IT$reg_id), cor(IT$nat_id,IT$eur_id), cor(IT$reg_id,IT$eur_id)), digits=2)
# ES
Table2[4,] <-round(c(mean(ES$pro_imm),sd(ES$pro_imm), mean(ES$pro_red),sd(ES$pro_red), cor(ES$pro_imm,ES$pro_red),
                     mean(ES$nat_id),sd(ES$nat_id), mean(ES$reg_id),sd(ES$reg_id), mean(ES$eur_id),sd(ES$eur_id),
                     cor(ES$nat_id,ES$reg_id), cor(ES$nat_id,ES$eur_id), cor(ES$reg_id,ES$eur_id)), digits=2)
# PL
Table2[5,] <-round(c(mean(PL$pro_imm),sd(PL$pro_imm), mean(PL$pro_red),sd(PL$pro_red), cor(PL$pro_imm,PL$pro_red),
                     mean(PL$nat_id),sd(PL$nat_id), mean(PL$reg_id),sd(PL$reg_id), mean(PL$eur_id),sd(PL$eur_id),
                     cor(PL$nat_id,PL$reg_id), cor(PL$nat_id,PL$eur_id), cor(PL$reg_id,PL$eur_id)), digits=2)
# RO
Table2[6,] <-round(c(mean(RO$pro_imm),sd(RO$pro_imm), mean(RO$pro_red),sd(RO$pro_red), cor(RO$pro_imm,RO$pro_red),
                     mean(RO$nat_id),sd(RO$nat_id), mean(RO$reg_id),sd(RO$reg_id), mean(RO$eur_id),sd(RO$eur_id),
                     cor(RO$nat_id,RO$reg_id), cor(RO$nat_id,RO$eur_id), cor(RO$reg_id,RO$eur_id)), digits=2)
# NL
Table2[7,]<-round(c(mean(NL$pro_imm),sd(NL$pro_imm), mean(NL$pro_red),sd(NL$pro_red), cor(NL$pro_imm,NL$pro_red),
                    mean(NL$nat_id),sd(NL$nat_id), mean(NL$reg_id),sd(NL$reg_id), mean(NL$eur_id),sd(NL$eur_id),
                    cor(NL$nat_id,NL$reg_id), cor(NL$nat_id,NL$eur_id), cor(NL$reg_id,NL$eur_id)), digits=2)
# SE
Table2[8,] <-round(c(mean(SE$pro_imm),sd(SE$pro_imm), mean(SE$pro_red),sd(SE$pro_red), cor(SE$pro_imm,SE$pro_red),
                     mean(SE$nat_id),sd(SE$nat_id), mean(SE$reg_id),sd(SE$reg_id), mean(SE$eur_id),sd(SE$eur_id),
                     cor(SE$nat_id,SE$reg_id), cor(SE$nat_id,SE$eur_id), cor(SE$reg_id,SE$eur_id)), digits=2)
# HU
Table2[9,] <-round(c(mean(HU$pro_imm),sd(HU$pro_imm), mean(HU$pro_red),sd(HU$pro_red), cor(HU$pro_imm,HU$pro_red),
                     mean(HU$nat_id),sd(HU$nat_id), mean(HU$reg_id),sd(HU$reg_id), mean(HU$eur_id), sd(HU$eur_id),
                     cor(HU$nat_id,HU$reg_id), cor(HU$nat_id,HU$eur_id), cor(HU$reg_id,HU$eur_id)), digits=2)
# AT
Table2[10,] <-round(c(mean(AT$pro_imm),  sd(AT$pro_imm), mean(AT$pro_red), sd(AT$pro_red), cor(AT$pro_imm,AT$pro_red),
                      mean(AT$nat_id),sd(AT$nat_id), mean(AT$reg_id),sd(AT$reg_id), mean(AT$eur_id),sd(AT$eur_id),
                      cor(AT$nat_id,AT$reg_id), cor(AT$nat_id,AT$eur_id), cor(AT$reg_id,AT$eur_id)), digits=2)
# BG
Table2[11,] <-round(c(mean(BG$pro_imm),sd(BG$pro_imm), mean(BG$pro_red),sd(BG$pro_red), cor(BG$pro_imm,BG$pro_red),
                      mean(BG$nat_id),sd(BG$nat_id), mean(BG$reg_id),sd(BG$reg_id), mean(BG$eur_id),sd(BG$eur_id),
                      cor(BG$nat_id,BG$reg_id), cor(BG$nat_id,BG$eur_id), cor(BG$reg_id,BG$eur_id)), digits=2)
# SK
Table2[12,] <-round(c(mean(SK$pro_imm),sd(SK$pro_imm), mean(SK$pro_red),sd(SK$pro_red), cor(SK$pro_imm,SK$pro_red),
                      mean(SK$nat_id),sd(SK$nat_id), mean(SK$reg_id),sd(SK$reg_id), mean(SK$eur_id),sd(SK$eur_id),
                      cor(SK$nat_id,SK$reg_id), cor(SK$nat_id,SK$eur_id), cor(SK$reg_id,SK$eur_id)), digits=2)
# LV
Table2[13,] <-round(c(mean(LV$pro_imm),sd(LV$pro_imm), mean(LV$pro_red),sd(LV$pro_red), cor(LV$pro_imm,LV$pro_red),
                      mean(LV$nat_id),sd(LV$nat_id), mean(LV$reg_id),sd(LV$reg_id), mean(LV$eur_id),sd(LV$eur_id),
                      cor(LV$nat_id,LV$reg_id), cor(LV$nat_id,LV$eur_id), cor(LV$reg_id,LV$eur_id)), digits=2)
# EE
Table2[14,] <-round(c(mean(EE$pro_imm),sd(EE$pro_imm), mean(EE$pro_red),sd(EE$pro_red), cor(EE$pro_imm,EE$pro_red),
                      mean(EE$nat_id),sd(EE$nat_id), mean(EE$reg_id),sd(EE$reg_id), mean(EE$eur_id),sd(EE$eur_id),
                      cor(EE$nat_id,EE$reg_id), cor(EE$nat_id,EE$eur_id), cor(EE$reg_id,EE$eur_id)), digits=2)

Table2


#-------------------------------------------------------------------
# Table 3: Descriptive Statistics on Control Variables
#--------------------------------------------------------------------
Table3 <- matrix(NA, 14, ncol = 13)
colnames(Table3) <- c("female","age","rural", "small","large","elem.","high s.","uni.", "econ.sat.",
                      "low","medium","high", "EU funds")
rownames(Table3) <- c("Germany", "France", "Italy", "Spain", "Poland", "Romania", "Netherlands", "Sweden", "Hungary", "Austria", "Belgium", "Slovakia", "Latvia", "Estonia")


# DE
Table3[1,] <- round(c(prop.table(table(DE$female))[2], mean(DE$age), 
                      prop.table(table(DE$pop))[1], prop.table(table(DE$pop))[2], prop.table(table(DE$pop))[3],
                      prop.table(table(DE$educ))[1], prop.table(table(DE$educ))[2], prop.table(table(DE$educ))[3],
                      prop.table(table(DE$eco))[2], prop.table(table(DE$income))[1], prop.table(table(DE$income))[2], prop.table(table(DE$income))[3],
                      prop.table(table(DE$eu_ben))[2]), digits=2)
# FR
Table3[2,] <- round(c(prop.table(table(FR$female))[2], mean(FR$age),
                      prop.table(table(FR$pop))[1], prop.table(table(FR$pop))[2], prop.table(table(FR$pop))[3],
                      prop.table(table(FR$educ))[1], prop.table(table(FR$educ))[2], prop.table(table(FR$educ))[3],
                      prop.table(table(FR$eco))[2], prop.table(table(FR$income))[1], prop.table(table(FR$income))[2], 
                      prop.table(table(FR$income))[3],prop.table(table(FR$eu_ben))[2]), digits=2)
# IT
Table3[3,] <- round(c(prop.table(table(IT$female))[2], mean(IT$age),
                      prop.table(table(IT$pop))[1], prop.table(table(IT$pop))[2], prop.table(table(IT$pop))[3],
                      prop.table(table(IT$educ))[1], prop.table(table(IT$educ))[2], prop.table(table(IT$educ))[3],
                      prop.table(table(IT$eco))[2], prop.table(table(IT$income))[1], prop.table(table(IT$income))[2], 
                      prop.table(table(IT$income))[3],prop.table(table(IT$eu_ben))[2]), digits=2)
# ES
Table3[4,] <- round(c(prop.table(table(ES$female))[2], mean(ES$age),
                      prop.table(table(ES$pop))[1], prop.table(table(ES$pop))[2], prop.table(table(ES$pop))[3],
                      prop.table(table(ES$educ))[1], prop.table(table(ES$educ))[2], prop.table(table(ES$educ))[3],
                      prop.table(table(ES$eco))[2], prop.table(table(ES$income))[1], prop.table(table(ES$income))[2], 
                      prop.table(table(ES$income))[3],prop.table(table(ES$eu_ben))[2]), digits=2)
# PL
Table3[5,] <- round(c(prop.table(table(PL$female))[2], mean(PL$age),
                      prop.table(table(PL$pop))[1], prop.table(table(PL$pop))[2], prop.table(table(PL$pop))[3],
                      prop.table(table(PL$educ))[1], prop.table(table(PL$educ))[2], prop.table(table(PL$educ))[3],
                      prop.table(table(PL$eco))[2], prop.table(table(PL$income))[1], prop.table(table(PL$income))[2], 
                      prop.table(table(PL$income))[3],prop.table(table(PL$eu_ben))[2]), digits=2)
# RO
Table3[6,] <- round(c(prop.table(table(RO$female))[2], mean(RO$age),
                      prop.table(table(RO$pop))[1], prop.table(table(RO$pop))[2], prop.table(table(RO$pop))[3],
                      prop.table(table(RO$educ))[1], prop.table(table(RO$educ))[2], prop.table(table(RO$educ))[3],
                      prop.table(table(RO$eco))[2], prop.table(table(RO$income))[1], prop.table(table(RO$income))[2], 
                      prop.table(table(RO$income))[3],prop.table(table(RO$eu_ben))[2]), digits=2)
# NL 
Table3[7,] <- round(c(prop.table(table(NL$female))[2], mean(NL$age),
                      prop.table(table(NL$pop))[1], prop.table(table(NL$pop))[2], prop.table(table(NL$pop))[3],
                      prop.table(table(NL$educ))[1], prop.table(table(NL$educ))[2], prop.table(table(NL$educ))[3],
                      prop.table(table(NL$eco))[2], prop.table(table(NL$income))[1], prop.table(table(NL$income))[2], 
                      prop.table(table(NL$income))[3],prop.table(table(NL$eu_ben))[2]), digits=2)
# SE
Table3[8,] <- round(c(prop.table(table(SE$female))[2], mean(SE$age),
                      prop.table(table(SE$pop))[1], prop.table(table(SE$pop))[2], prop.table(table(SE$pop))[3],
                      prop.table(table(SE$educ))[1], prop.table(table(SE$educ))[2], prop.table(table(SE$educ))[3],
                      prop.table(table(SE$eco))[2], prop.table(table(SE$income))[1], prop.table(table(SE$income))[2], 
                      prop.table(table(SE$income))[3],prop.table(table(SE$eu_ben))[2]), digits=2)
# HU
Table3[9,] <- round(c(prop.table(table(HU$female))[2], mean(HU$age),
                      prop.table(table(HU$pop))[1], prop.table(table(HU$pop))[2], prop.table(table(HU$pop))[3],
                      prop.table(table(HU$educ))[1], prop.table(table(HU$educ))[2], prop.table(table(HU$educ))[3],
                      prop.table(table(HU$eco))[2], prop.table(table(HU$income))[1], prop.table(table(HU$income))[2], 
                      prop.table(table(HU$income))[3],prop.table(table(HU$eu_ben))[2]), digits=2)
# AT
Table3[10,] <- round(c(prop.table(table(AT$female))[2], mean(AT$age),
                       prop.table(table(AT$pop))[1], prop.table(table(AT$pop))[2], prop.table(table(AT$pop))[3],
                       prop.table(table(AT$educ))[1], prop.table(table(AT$educ))[2], prop.table(table(AT$educ))[3],
                       prop.table(table(AT$eco))[2], prop.table(table(AT$income))[1], prop.table(table(AT$income))[2], 
                       prop.table(table(AT$income))[3],prop.table(table(AT$eu_ben))[2]), digits=2)
# BG
Table3[11,] <- round(c(prop.table(table(BG$female))[2], mean(BG$age),
                       prop.table(table(BG$pop))[1], prop.table(table(BG$pop))[2], prop.table(table(BG$pop))[3],
                       prop.table(table(BG$educ))[1], prop.table(table(BG$educ))[2], prop.table(table(BG$educ))[3],
                       prop.table(table(BG$eco))[2], prop.table(table(BG$income))[1], prop.table(table(BG$income))[2], 
                       prop.table(table(BG$income))[3],prop.table(table(BG$eu_ben))[2]), digits=2)
# SK 
Table3[12,] <- round(c(prop.table(table(SK$female))[2], mean(SK$age),
                       prop.table(table(SK$pop))[1], prop.table(table(SK$pop))[2], prop.table(table(SK$pop))[3],
                       prop.table(table(SK$educ))[1], prop.table(table(SK$educ))[2], prop.table(table(SK$educ))[3],
                       prop.table(table(SK$eco))[2], prop.table(table(SK$income))[1], prop.table(table(SK$income))[2], 
                       prop.table(table(SK$income))[3],prop.table(table(SK$eu_ben))[2]), digits=2)
# LV 
Table3[13,] <- round(c(prop.table(table(LV$female))[2], mean(LV$age),
                       prop.table(table(LV$pop))[1], prop.table(table(LV$pop))[2], prop.table(table(LV$pop))[3],
                       prop.table(table(LV$educ))[1], prop.table(table(LV$educ))[2], prop.table(table(LV$educ))[3],
                       prop.table(table(LV$eco))[2], prop.table(table(LV$income))[1], prop.table(table(LV$income))[2], 
                       prop.table(table(LV$income))[3],prop.table(table(LV$eu_ben))[2]), digits=2)
# EE
Table3[14,] <- round(c(prop.table(table(EE$female))[2], mean(EE$age),
                       prop.table(table(EE$pop))[1], prop.table(table(EE$pop))[2], prop.table(table(EE$pop))[3],
                       prop.table(table(EE$educ))[1], prop.table(table(EE$educ))[2], prop.table(table(EE$educ))[3],
                       prop.table(table(EE$eco))[2], prop.table(table(EE$income))[1], prop.table(table(EE$income))[2], 
                       prop.table(table(EE$income))[3],prop.table(table(EE$eu_ben))[2]), digits=2)

Table3


#===================================================================================
# Estimations (for Table 4 and Table 5)
#===================================================================================

#---------------------------
# Germany (DE)
#---------------------------
DE_1 <- glm(eu_support ~ female + aged + pop + educ + eco + income + eu_ben, data=DE, family="binomial")
DE_2 <- glm(eu_support ~ pro_imms + pro_reds + female + aged + pop + educ + eco + income + eu_ben, data=DE, family="binomial")
DE_3 <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, data=DE, family="binomial")
#---------------------------
# France (FR)
#---------------------------
FR_1 <- glm(eu_support ~ female + aged + pop + educ + eco + income + eu_ben, data=FR, family="binomial")
FR_2 <- glm(eu_support ~ pro_imms + pro_reds + female + aged + pop + educ + eco + income + eu_ben, data=FR, family="binomial")
FR_3 <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, data=FR, family="binomial")
#---------------------------
# Italy (IT)
#---------------------------
IT_1 <- glm(eu_support ~ female + aged + pop + educ + eco + income + eu_ben, data=IT, family="binomial")
IT_2 <- glm(eu_support ~ pro_imms + pro_reds + female + aged + pop + educ + eco + income + eu_ben, data=IT, family="binomial")
IT_3 <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, data=IT, family="binomial")
#---------------------------
# Spain (ES)
#---------------------------
ES_1 <- glm(eu_support ~ female + aged + pop + educ + eco + income + eu_ben, data=ES, family="binomial")
ES_2 <- glm(eu_support ~ pro_imms + pro_reds + female + aged + pop + educ + eco + income + eu_ben, data=ES, family="binomial")
ES_3 <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, data=ES, family="binomial")
#---------------------------
# Poland (PL)
#---------------------------
PL_1 <- glm(eu_support ~ female + aged + pop + educ + eco + income + eu_ben, data=PL, family="binomial")
PL_2 <- glm(eu_support ~ pro_imms + pro_reds + female + aged + pop + educ + eco + income + eu_ben, data=PL, family="binomial")
PL_3 <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, data=PL, family="binomial")
#---------------------------
# Romania (RO)
#---------------------------
RO_1 <- glm(eu_support ~ female + aged + pop + educ + eco + income + eu_ben, data=RO, family="binomial")
RO_2 <- glm(eu_support ~ pro_imms + pro_reds + female + aged + pop + educ + eco + income + eu_ben, data=RO, family="binomial")
RO_3 <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, data=RO, family="binomial")
#---------------------------
# Netherlands (NL)
#---------------------------
NL_1 <- glm(eu_support ~ female + aged + pop + educ + eco + income + eu_ben, data=NL, family="binomial")
NL_2 <- glm(eu_support ~ pro_imms + pro_reds + female + aged + pop + educ + eco + income + eu_ben, data=NL, family="binomial")
NL_3 <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, data=NL, family="binomial")
#---------------------------
# Sweden (SE)
#---------------------------
SE_1 <- glm(eu_support ~ female + aged + pop + educ + eco + income + eu_ben, data=SE, family="binomial")
SE_2 <- glm(eu_support ~ pro_imms + pro_reds + female + aged + pop + educ + eco + income + eu_ben, data=SE, family="binomial")
SE_3 <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, data=SE, family="binomial")
#---------------------------
# Hungary (HU)
#---------------------------
HU_1 <- glm(eu_support ~ female + aged + pop + educ + eco + income + eu_ben, data=HU, family="binomial")
HU_2 <- glm(eu_support ~ pro_imms + pro_reds + female + aged + pop + educ + eco + income + eu_ben, data=HU, family="binomial")
HU_3 <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, data=HU, family="binomial")
#---------------------------
# Austria (AT)
#---------------------------
AT_1 <- glm(eu_support ~ female + aged + pop + educ + eco + income + eu_ben, data=AT, family="binomial")
AT_2 <- glm(eu_support ~ pro_imms + pro_reds + female + aged + pop + educ + eco + income + eu_ben, data=AT, family="binomial")
AT_3 <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, data=AT, family="binomial")
#---------------------------
# Bulgaria (BG)
#---------------------------
BG_1 <- glm(eu_support ~ female + aged + pop + educ + eco + income + eu_ben, data=BG, family="binomial")
BG_2 <- glm(eu_support ~ pro_imms + pro_reds + female + aged + pop + educ + eco + income + eu_ben, data=BG, family="binomial")
BG_3 <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, data=BG, family="binomial")
#---------------------------
# Slovakia (SK)
#---------------------------
SK_1 <- glm(eu_support ~ female + aged + pop + educ + eco + income + eu_ben, data=SK, family="binomial")
SK_2 <- glm(eu_support ~ pro_imms + pro_reds + female + aged + pop + educ + eco + income + eu_ben, data=SK, family="binomial")
SK_3 <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, data=SK, family="binomial")
#---------------------------
# Latvia (LV)
#---------------------------
LV_1 <- glm(eu_support ~ female + aged + pop + educ + eco + income + eu_ben, data=LV, family="binomial")
LV_2 <- glm(eu_support ~ pro_imms + pro_reds + female + aged + pop + educ + eco + income + eu_ben, data=LV, family="binomial")
LV_3 <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, data=LV, family="binomial")
#---------------------------
# Estonia (EE)
#---------------------------
EE_1 <- glm(eu_support ~ female + aged + pop + educ + eco + income + eu_ben, data=EE, family="binomial")
EE_2 <- glm(eu_support ~ pro_imms + pro_reds + female + aged + pop + educ + eco + income + eu_ben, data=EE, family="binomial")
EE_3 <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, data=EE, family="binomial")
#---------------------------
# Pooled
#---------------------------
FS_P <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, data=FS, family="binomial")


#-------------------------------------------------------------------
# Table 4: Explaining Support for EU Membership: Model Comparisons
#--------------------------------------------------------------------
# AIC
Table4 <- matrix(NA, 14, ncol = 4)
colnames(Table4) <- c("Obs.","M1: Controls z_i", "M2: M1 + x_i + y_i", "M3: M2 + (n-e) + (n-r)")
rownames(Table4) <- c("Germany", "France", "Italy", "Spain", "Poland", "Romania", "Netherlands",
                      "Sweden", "Hungary","Austria",  "Bulgaria", "Slovakia", "Latvia", "Estonia")
Table4[1,]   <- c(nrow(DE),round(c(AIC(DE_1), AIC(DE_2), AIC(DE_3)), digits=2))
Table4[2,]   <- c(nrow(FR),round(c(AIC(FR_1), AIC(FR_2), AIC(FR_3)), digits=2))
Table4[3,]   <- c(nrow(IT),round(c(AIC(IT_1), AIC(IT_2), AIC(IT_3)), digits=2))
Table4[4,]   <- c(nrow(ES),round(c(AIC(ES_1), AIC(ES_2), AIC(ES_3)), digits=2))
Table4[5,]   <- c(nrow(PL),round(c(AIC(PL_1), AIC(PL_2), AIC(PL_3)), digits=2))
Table4[6,]   <- c(nrow(RO),round(c(AIC(RO_1), AIC(RO_2), AIC(RO_3)), digits=2))
Table4[7,]   <- c(nrow(NL),round(c(AIC(NL_1), AIC(NL_2), AIC(NL_3)), digits=2))
Table4[8,]   <- c(nrow(SE),round(c(AIC(SE_1), AIC(SE_2), AIC(SE_3)), digits=2))
Table4[9,]   <- c(nrow(HU),round(c(AIC(HU_1), AIC(HU_2), AIC(HU_3)), digits=2))
Table4[10,]  <- c(nrow(AT),round(c(AIC(AT_1), AIC(AT_2), AIC(AT_3)), digits=2))
Table4[11,]  <- c(nrow(BG),round(c(AIC(BG_1), AIC(BG_2), AIC(BG_3)), digits=2))
Table4[12,]  <- c(nrow(SK),round(c(AIC(SK_1), AIC(SK_2), AIC(SK_3)), digits=2))
Table4[13,]  <- c(nrow(LV),round(c(AIC(LV_1), AIC(LV_2), AIC(LV_3)), digits=2))
Table4[14,]  <- c(nrow(EE),round(c(AIC(EE_1), AIC(EE_2), AIC(EE_3)), digits=2))

Table4


#---------------------------------------------------------------------------------------
# Table 5: Effects of Ideological Variables and Territorial Identity Differentials (M3)
#----------------------------------------------------------------------------------------
Table5 <- matrix(NA, 15, ncol = 10)
colnames(Table5) <- c("beta_1", "", "beta_2", "", "beta_3", "", "beta_4", "", "Log-Lik.", "R^2")
rownames(Table5) <- c("Germany", "France", "Italy", "Spain", "Poland", "Romania", "Netherlands",
                      "Sweden", "Hungary","Austria",  "Bulgaria", "Slovakia", "Latvia", "Estonia","Pooled")

# DE
Table5[1,] <-round(c(summary(DE_3)$coefficients[2,1], summary(DE_3)$coefficients[2,2],
                     summary(DE_3)$coefficients[3,1], summary(DE_3)$coefficients[3,2],
                     summary(DE_3)$coefficients[4,1], summary(DE_3)$coefficients[4,2],
                     summary(DE_3)$coefficients[5,1], summary(DE_3)$coefficients[5,2],
                     logLik(DE_3), with(summary(DE_3), 1-deviance/null.deviance)),digits=2)
# FR
Table5[2,] <-round(c(summary(FR_3)$coefficients[2,1], summary(FR_3)$coefficients[2,2],
                     summary(FR_3)$coefficients[3,1], summary(FR_3)$coefficients[3,2],
                     summary(FR_3)$coefficients[4,1], summary(FR_3)$coefficients[4,2],
                     summary(FR_3)$coefficients[5,1], summary(FR_3)$coefficients[5,2],
                     logLik(FR_3), with(summary(FR_3), 1-deviance/null.deviance)),digits=2)
# IT
Table5[3,] <-round(c(summary(IT_3)$coefficients[2,1], summary(IT_3)$coefficients[2,2],
                     summary(IT_3)$coefficients[3,1], summary(IT_3)$coefficients[3,2],
                     summary(IT_3)$coefficients[4,1], summary(IT_3)$coefficients[4,2],
                     summary(IT_3)$coefficients[5,1], summary(IT_3)$coefficients[5,2],
                     logLik(IT_3), with(summary(IT_3), 1-deviance/null.deviance)),digits=2)
# ES
Table5[4,] <-round(c(summary(ES_3)$coefficients[2,1], summary(ES_3)$coefficients[2,2],
                     summary(ES_3)$coefficients[3,1], summary(ES_3)$coefficients[3,2],
                     summary(ES_3)$coefficients[4,1], summary(ES_3)$coefficients[4,2],
                     summary(ES_3)$coefficients[5,1], summary(ES_3)$coefficients[5,2],
                     logLik(ES_3), with(summary(ES_3), 1-deviance/null.deviance)),digits=2)
# PL
Table5[5,] <-round(c(summary(PL_3)$coefficients[2,1], summary(PL_3)$coefficients[2,2],
                     summary(PL_3)$coefficients[3,1], summary(PL_3)$coefficients[3,2],
                     summary(PL_3)$coefficients[4,1], summary(PL_3)$coefficients[4,2],
                     summary(PL_3)$coefficients[5,1], summary(PL_3)$coefficients[5,2],
                     logLik(PL_3), with(summary(PL_3), 1-deviance/null.deviance)),digits=2)
# RO
Table5[6,] <-round(c(summary(RO_3)$coefficients[2,1], summary(RO_3)$coefficients[2,2],
                     summary(RO_3)$coefficients[3,1], summary(RO_3)$coefficients[3,2],
                     summary(RO_3)$coefficients[4,1], summary(RO_3)$coefficients[4,2],
                     summary(RO_3)$coefficients[5,1], summary(RO_3)$coefficients[5,2],
                     logLik(RO_3), with(summary(RO_3), 1-deviance/null.deviance)),digits=2)
# NL
Table5[7,] <-round(c(summary(NL_3)$coefficients[2,1], summary(NL_3)$coefficients[2,2],
                     summary(NL_3)$coefficients[3,1], summary(NL_3)$coefficients[3,2],
                     summary(NL_3)$coefficients[4,1], summary(NL_3)$coefficients[4,2],
                     summary(NL_3)$coefficients[5,1], summary(NL_3)$coefficients[5,2],
                     logLik(NL_3), with(summary(NL_3), 1-deviance/null.deviance)),digits=2)
# SE
Table5[8,] <-round(c(summary(SE_3)$coefficients[2,1], summary(SE_3)$coefficients[2,2],
                     summary(SE_3)$coefficients[3,1], summary(SE_3)$coefficients[3,2],
                     summary(SE_3)$coefficients[4,1], summary(SE_3)$coefficients[4,2],
                     summary(SE_3)$coefficients[5,1], summary(SE_3)$coefficients[5,2],
                     logLik(SE_3), with(summary(SE_3), 1-deviance/null.deviance)),digits=2)
# HU
Table5[9,] <-round(c(summary(HU_3)$coefficients[2,1], summary(HU_3)$coefficients[2,2],
                     summary(HU_3)$coefficients[3,1], summary(HU_3)$coefficients[3,2],
                     summary(HU_3)$coefficients[4,1], summary(HU_3)$coefficients[4,2],
                     summary(HU_3)$coefficients[5,1], summary(HU_3)$coefficients[5,2],
                     logLik(HU_3), with(summary(HU_3), 1-deviance/null.deviance)),digits=2)
# AT
Table5[10,] <-round(c(summary(AT_3)$coefficients[2,1], summary(AT_3)$coefficients[2,2],
                      summary(AT_3)$coefficients[3,1], summary(AT_3)$coefficients[3,2],
                      summary(AT_3)$coefficients[4,1], summary(AT_3)$coefficients[4,2],
                      summary(AT_3)$coefficients[5,1], summary(AT_3)$coefficients[5,2],
                      logLik(AT_3), with(summary(AT_3), 1-deviance/null.deviance)),digits=2)
# BG
Table5[11,] <-round(c(summary(BG_3)$coefficients[2,1], summary(BG_3)$coefficients[2,2],
                      summary(BG_3)$coefficients[3,1], summary(BG_3)$coefficients[3,2],
                      summary(BG_3)$coefficients[4,1], summary(BG_3)$coefficients[4,2],
                      summary(BG_3)$coefficients[5,1], summary(BG_3)$coefficients[5,2],
                      logLik(BG_3), with(summary(BG_3), 1-deviance/null.deviance)),digits=2)
# SK
Table5[12,] <-round(c(summary(SK_3)$coefficients[2,1], summary(SK_3)$coefficients[2,2],
                      summary(SK_3)$coefficients[3,1], summary(SK_3)$coefficients[3,2],
                      summary(SK_3)$coefficients[4,1], summary(SK_3)$coefficients[4,2],
                      summary(SK_3)$coefficients[5,1], summary(SK_3)$coefficients[5,2],
                      logLik(SK_3), with(summary(SK_3), 1-deviance/null.deviance)),digits=2)
# LV
Table5[13,] <-round(c(summary(LV_3)$coefficients[2,1], summary(LV_3)$coefficients[2,2],
                      summary(LV_3)$coefficients[3,1], summary(LV_3)$coefficients[3,2],
                      summary(LV_3)$coefficients[4,1], summary(LV_3)$coefficients[4,2],
                      summary(LV_3)$coefficients[5,1], summary(LV_3)$coefficients[5,2],
                      logLik(LV_3), with(summary(LV_3), 1-deviance/null.deviance)),digits=2)
# EE
Table5[14,] <-round(c(summary(EE_3)$coefficients[2,1], summary(EE_3)$coefficients[2,2],
                      summary(EE_3)$coefficients[3,1], summary(EE_3)$coefficients[3,2],
                      summary(EE_3)$coefficients[4,1], summary(EE_3)$coefficients[4,2],
                      summary(EE_3)$coefficients[5,1], summary(EE_3)$coefficients[5,2],
                      logLik(EE_3), with(summary(EE_3), 1-deviance/null.deviance)),digits=2)

# POOLED
Table5[15,] <-round(c(summary(FS_P)$coefficients[2,1], summary(FS_P)$coefficients[2,2],
                      summary(FS_P)$coefficients[3,1], summary(FS_P)$coefficients[3,2],
                      summary(FS_P)$coefficients[4,1], summary(FS_P)$coefficients[4,2],
                      summary(FS_P)$coefficients[5,1], summary(FS_P)$coefficients[5,2],
                      logLik(FS_P), with(summary(FS_P), 1-deviance/null.deviance)),digits=2)

Table5


#---------------------------------------------------------------------------------------
# Figure 3: Impact of Ideological Considerations and Territorial Identity Differentials
#---------------------------------------------------------------------------------------

#=======================================================
# Immigration
#=======================================================

pdf(file="Figure_3a.pdf", width=10, height=9)


country <- rep(c("Germany","France","Italy","Spain","Poland","Romania","Netherlands","Sweden","Hungary","Austria",
                 "Bulgaria","Slovakia","Latvia","Estonia","Pooled"))  
coef <- c(round(summary(DE_3)$coefficients[2,1], digits=2),
          round(summary(FR_3)$coefficients[2,1], digits=2),
          round(summary(IT_3)$coefficients[2,1], digits=2),
          round(summary(ES_3)$coefficients[2,1], digits=2),
          round(summary(PL_3)$coefficients[2,1], digits=2),
          round(summary(RO_3)$coefficients[2,1], digits=2),
          round(summary(NL_3)$coefficients[2,1], digits=2),
          round(summary(SE_3)$coefficients[2,1], digits=2),
          round(summary(HU_3)$coefficients[2,1], digits=2),
          round(summary(AT_3)$coefficients[2,1], digits=2),
          round(summary(BG_3)$coefficients[2,1], digits=2),
          round(summary(SK_3)$coefficients[2,1], digits=2),
          round(summary(LV_3)$coefficients[2,1], digits=2),
          round(summary(EE_3)$coefficients[2,1], digits=2),
          round(summary(FS_P)$coefficients[2,1], digits=2))

lower <- c(confintvglm(DE_3)[2,1],
           confintvglm(FR_3)[2,1],
           confintvglm(IT_3)[2,1],
           confintvglm(ES_3)[2,1],
           confintvglm(PL_3)[2,1],
           confintvglm(RO_3)[2,1],
           confintvglm(NL_3)[2,1],
           confintvglm(SE_3)[2,1],
           confintvglm(HU_3)[2,1],
           confintvglm(AT_3)[2,1],
           confintvglm(BG_3)[2,1],
           confintvglm(SK_3)[2,1],
           confintvglm(LV_3)[2,1],
           confintvglm(EE_3)[2,1],
           confintvglm(FS_P)[2,1])

upper <- c(confintvglm(DE_3)[2,2],
           confintvglm(FR_3)[2,2],
           confintvglm(IT_3)[2,2],
           confintvglm(ES_3)[2,2],
           confintvglm(PL_3)[2,2],
           confintvglm(RO_3)[2,2],
           confintvglm(NL_3)[2,2],
           confintvglm(SE_3)[2,2],
           confintvglm(HU_3)[2,2],
           confintvglm(AT_3)[2,2],
           confintvglm(BG_3)[2,2],
           confintvglm(SK_3)[2,2],
           confintvglm(LV_3)[2,2],
           confintvglm(EE_3)[2,2],
           confintvglm(FS_P)[2,2])


df <- data.frame(country = country,
                 coefficients = coef,
                 IC_lower = lower,
                 IC_upper = upper)

df$country <- factor(df$country,
                     levels = c("Pooled","Germany","France","Italy","Spain","Poland","Romania","Netherlands",
                                "Sweden","Hungary","Austria","Bulgaria","Slovakia","Latvia","Estonia"))

shapes <- c(16, 17, 15, 3, 4, 18, 8, 7, 9, 10, 12, 13, 14, 0, 1)
colors <- c("#000000", "#1f78b4", "#e41a1c", "#4daf4a", "#ff7f00",
            "#984ea3", "#377eb8", "#a65628", "#f781bf", "#999999",
            "#66c2a5", "#fc8d62", "#8da0cb", "#e78ac3", "#a6d854")

ggplot(df, aes(x = coefficients,
               y = country,
               xmin = IC_lower,
               xmax = IC_upper,
               color = country,
               shape = country)) + 
  geom_vline(xintercept = seq(-2, 2, by = 0.5), linewidth = 1.0, color = "grey")+
  geom_vline(xintercept = 0, linewidth = 1.2, linetype = "dashed", color = "black")+
  geom_pointrange(position = position_dodge(width = 0.7),
                  size = 1.5) +
  scale_color_manual(values = colors) +
  scale_shape_manual(values = shapes) +
  coord_flip() +
  theme_minimal(base_size = 13) +
  theme(axis.text.y = element_text(size = 20),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 20),
        legend.position = "bottom",
        legend.title = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.major.x = element_line(linewidth = 1.5)
  ) +
  
  labs(
    x = NULL,
    y = NULL
  )+
  theme(legend.position = "none")

dev.off()

#=======================================================
# Redistribution
#=======================================================

pdf(file="Figure_3b.pdf", width=10, height=9)

country <- rep(c("Germany","France","Italy","Spain","Poland","Romania","Netherlands","Sweden","Hungary","Austria","Bulgaria","Slovakia","Latvia","Estonia","Pooled"))  
coef <- c(round(summary(DE_3)$coefficients[3,1], digits=2),
          round(summary(FR_3)$coefficients[3,1], digits=2),
          round(summary(IT_3)$coefficients[3,1], digits=2),
          round(summary(ES_3)$coefficients[3,1], digits=2),
          round(summary(PL_3)$coefficients[3,1], digits=2),
          round(summary(RO_3)$coefficients[3,1], digits=2),
          round(summary(NL_3)$coefficients[3,1], digits=2),
          round(summary(SE_3)$coefficients[3,1], digits=2),
          round(summary(HU_3)$coefficients[3,1], digits=2),
          round(summary(AT_3)$coefficients[3,1], digits=2),
          round(summary(BG_3)$coefficients[3,1], digits=2),
          round(summary(SK_3)$coefficients[3,1], digits=2),
          round(summary(LV_3)$coefficients[3,1], digits=2),
          round(summary(EE_3)$coefficients[3,1], digits=2),
          round(summary(FS_P)$coefficients[3,1], digits=2))

lower <- c(confintvglm(DE_3)[3,1],
           confintvglm(FR_3)[3,1],
           confintvglm(IT_3)[3,1],
           confintvglm(ES_3)[3,1],
           confintvglm(PL_3)[3,1],
           confintvglm(RO_3)[3,1],
           confintvglm(NL_3)[3,1],
           confintvglm(SE_3)[3,1],
           confintvglm(HU_3)[3,1],
           confintvglm(AT_3)[3,1],
           confintvglm(BG_3)[3,1],
           confintvglm(SK_3)[3,1],
           confintvglm(LV_3)[3,1],
           confintvglm(EE_3)[3,1],
           confintvglm(FS_P)[3,1])

upper <- c(confintvglm(DE_3)[3,2],
           confintvglm(FR_3)[3,2],
           confintvglm(IT_3)[3,2],
           confintvglm(ES_3)[3,2],
           confintvglm(PL_3)[3,2],
           confintvglm(RO_3)[3,2],
           confintvglm(NL_3)[3,2],
           confintvglm(SE_3)[3,2],
           confintvglm(HU_3)[3,2],
           confintvglm(AT_3)[3,2],
           confintvglm(BG_3)[3,2],
           confintvglm(SK_3)[3,2],
           confintvglm(LV_3)[3,2],
           confintvglm(EE_3)[3,2],
           confintvglm(FS_P)[3,2])


df <- data.frame(country = country,
                 coefficients = coef,
                 IC_lower = lower,
                 IC_upper = upper)

df$country <- factor(df$country,
                     levels = c("Pooled","Germany","France","Italy","Spain","Poland","Romania","Netherlands",
                                "Sweden","Hungary","Austria","Bulgaria","Slovakia","Latvia","Estonia"))

shapes <- c(16, 17, 15, 3, 4, 18, 8, 7, 9, 10, 12, 13, 14, 0, 1)
colors <- c("#000000", "#1f78b4", "#e41a1c", "#4daf4a", "#ff7f00",
            "#984ea3", "#377eb8", "#a65628", "#f781bf", "#999999",
            "#66c2a5", "#fc8d62", "#8da0cb", "#e78ac3", "#a6d854")

ggplot(df, aes(x = coefficients,
               y = country,
               xmin = IC_lower,
               xmax = IC_upper,
               color = country,
               shape = country)) + 
  geom_vline(xintercept = seq(-2, 2, by = 0.5), linewidth = 1.0, color = "grey")+
  geom_vline(xintercept = 0, linewidth = 1.2, linetype = "dashed", color = "black")+
  geom_pointrange(position = position_dodge(width = 0.7),
                  size = 1.5) +
  scale_color_manual(values = colors) +
  scale_shape_manual(values = shapes) +
  coord_flip() +
  theme_minimal(base_size = 13) +
  theme(axis.text.y = element_text(size = 20),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 20),
        legend.position = "bottom",
        legend.title = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.major.x = element_line(linewidth = 1.5)
  ) +
  
  labs(
    x = NULL,
    y = NULL
  )+
  theme(legend.position = "none")


dev.off()

#=======================================================
# National-European
#=======================================================

pdf(file="Figure_3c.pdf", width=10, height=9)

country <- rep(c("Germany","France","Italy","Spain","Poland","Romania","Netherlands","Sweden","Hungary","Austria","Bulgaria","Slovakia","Latvia","Estonia","Pooled"))  
coef <- c(round(summary(DE_3)$coefficients[4,1], digits=2),
          round(summary(FR_3)$coefficients[4,1], digits=2),
          round(summary(IT_3)$coefficients[4,1], digits=2),
          round(summary(ES_3)$coefficients[4,1], digits=2),
          round(summary(PL_3)$coefficients[4,1], digits=2),
          round(summary(RO_3)$coefficients[4,1], digits=2),
          round(summary(NL_3)$coefficients[4,1], digits=2),
          round(summary(SE_3)$coefficients[4,1], digits=2),
          round(summary(HU_3)$coefficients[4,1], digits=2),
          round(summary(AT_3)$coefficients[4,1], digits=2),
          round(summary(BG_3)$coefficients[4,1], digits=2),
          round(summary(SK_3)$coefficients[4,1], digits=2),
          round(summary(LV_3)$coefficients[4,1], digits=2),
          round(summary(EE_3)$coefficients[4,1], digits=2),
          round(summary(FS_P)$coefficients[4,1], digits=2))

lower <- c(confintvglm(DE_3)[4,1],
           confintvglm(FR_3)[4,1],
           confintvglm(IT_3)[4,1],
           confintvglm(ES_3)[4,1],
           confintvglm(PL_3)[4,1],
           confintvglm(RO_3)[4,1],
           confintvglm(NL_3)[4,1],
           confintvglm(SE_3)[4,1],
           confintvglm(HU_3)[4,1],
           confintvglm(AT_3)[4,1],
           confintvglm(BG_3)[4,1],
           confintvglm(SK_3)[4,1],
           confintvglm(LV_3)[4,1],
           confintvglm(EE_3)[4,1],
           confintvglm(FS_P)[4,1])

upper <- c(confintvglm(DE_3)[4,2],
           confintvglm(FR_3)[4,2],
           confintvglm(IT_3)[4,2],
           confintvglm(ES_3)[4,2],
           confintvglm(PL_3)[4,2],
           confintvglm(RO_3)[4,2],
           confintvglm(NL_3)[4,2],
           confintvglm(SE_3)[4,2],
           confintvglm(HU_3)[4,2],
           confintvglm(AT_3)[4,2],
           confintvglm(BG_3)[4,2],
           confintvglm(SK_3)[4,2],
           confintvglm(LV_3)[4,2],
           confintvglm(EE_3)[4,2],
           confintvglm(FS_P)[4,2])


df <- data.frame(country = country,
                 coefficients = coef,
                 IC_lower = lower,
                 IC_upper = upper)

df$country <- factor(df$country,
                     levels = c("Pooled","Germany","France","Italy","Spain","Poland","Romania","Netherlands","Sweden",
                                "Hungary","Austria","Bulgaria","Slovakia","Latvia","Estonia"))

shapes <- c(16, 17, 15, 3, 4, 18, 8, 7, 9, 10, 12, 13, 14, 0, 1)
colors <- c("#000000", "#1f78b4", "#e41a1c", "#4daf4a", "#ff7f00",
            "#984ea3", "#377eb8", "#a65628", "#f781bf", "#999999",
            "#66c2a5", "#fc8d62", "#8da0cb", "#e78ac3", "#a6d854")

ggplot(df, aes(x = coefficients,
               y = country,
               xmin = IC_lower,
               xmax = IC_upper,
               color = country,
               shape = country)) + 
  geom_vline(xintercept = seq(-2, 2, by = 0.5), linewidth = 1.0, color = "grey")+
  geom_vline(xintercept = 0, linewidth = 1.2, linetype = "dashed", color = "black")+
  geom_pointrange(position = position_dodge(width = 0.7),
                  size = 1.5) +
  scale_color_manual(values = colors) +
  scale_shape_manual(values = shapes) +
  coord_flip() +
  theme_minimal(base_size = 13) +
  theme(axis.text.y = element_text(size = 20),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 20),
        legend.position = "bottom",
        legend.title = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.major.x = element_line(linewidth = 1.5)
  ) +
  
  labs(
    x = NULL,
    y = NULL
  )+
  theme(legend.position = "none")

dev.off()

#=======================================================
# National-Regional
#=======================================================

pdf(file="Figure_3d.pdf", width=10, height=9)

country <- rep(c("Germany","France","Italy","Spain","Poland","Romania","Netherlands","Sweden","Hungary","Austria","Bulgaria","Slovakia","Latvia","Estonia","Pooled"))  
coef <- c(round(summary(DE_3)$coefficients[5,1], digits=2),
          round(summary(FR_3)$coefficients[5,1], digits=2),
          round(summary(IT_3)$coefficients[5,1], digits=2),
          round(summary(ES_3)$coefficients[5,1], digits=2),
          round(summary(PL_3)$coefficients[5,1], digits=2),
          round(summary(RO_3)$coefficients[5,1], digits=2),
          round(summary(NL_3)$coefficients[5,1], digits=2),
          round(summary(SE_3)$coefficients[5,1], digits=2),
          round(summary(HU_3)$coefficients[5,1], digits=2),
          round(summary(AT_3)$coefficients[5,1], digits=2),
          round(summary(BG_3)$coefficients[5,1], digits=2),
          round(summary(SK_3)$coefficients[5,1], digits=2),
          round(summary(LV_3)$coefficients[5,1], digits=2),
          round(summary(EE_3)$coefficients[5,1], digits=2),
          round(summary(FS_P)$coefficients[5,1], digits=2))

lower <- c(confintvglm(DE_3)[5,1],
           confintvglm(FR_3)[5,1],
           confintvglm(IT_3)[5,1],
           confintvglm(ES_3)[5,1],
           confintvglm(PL_3)[5,1],
           confintvglm(RO_3)[5,1],
           confintvglm(NL_3)[5,1],
           confintvglm(SE_3)[5,1],
           confintvglm(HU_3)[5,1],
           confintvglm(AT_3)[5,1],
           confintvglm(BG_3)[5,1],
           confintvglm(SK_3)[5,1],
           confintvglm(LV_3)[5,1],
           confintvglm(EE_3)[5,1],
           confintvglm(FS_P)[5,1])

upper <- c(confintvglm(DE_3)[5,2],
           confintvglm(FR_3)[5,2],
           confintvglm(IT_3)[5,2],
           confintvglm(ES_3)[5,2],
           confintvglm(PL_3)[5,2],
           confintvglm(RO_3)[5,2],
           confintvglm(NL_3)[5,2],
           confintvglm(SE_3)[5,2],
           confintvglm(HU_3)[5,2],
           confintvglm(AT_3)[5,2],
           confintvglm(BG_3)[5,2],
           confintvglm(SK_3)[5,2],
           confintvglm(LV_3)[5,2],
           confintvglm(EE_3)[5,2],
           confintvglm(FS_P)[5,2])

df <- data.frame(country = country,
                 coefficients = coef,
                 IC_lower = lower,
                 IC_upper = upper)

df$country <- factor(df$country,
                     levels = c("Pooled","Germany","France","Italy","Spain","Poland","Romania","Netherlands","Sweden","Hungary",
                                "Austria","Bulgaria","Slovakia","Latvia","Estonia"))

shapes <- c(16, 17, 15, 3, 4, 18, 8, 7, 9, 10, 12, 13, 14, 0, 1)
colors <- c("#000000", "#1f78b4", "#e41a1c", "#4daf4a", "#ff7f00",
            "#984ea3", "#377eb8", "#a65628", "#f781bf", "#999999",
            "#66c2a5", "#fc8d62", "#8da0cb", "#e78ac3", "#a6d854")

ggplot(df, aes(x = coefficients,
               y = country,
               xmin = IC_lower,
               xmax = IC_upper,
               color = country,
               shape = country)) + 
  geom_vline(xintercept = seq(-2, 2, by = 0.5), linewidth = 1.0, color = "grey")+
  geom_vline(xintercept = 0, linewidth = 1.2, linetype = "dashed", color = "black")+
  geom_pointrange(position = position_dodge(width = 0.7),
                  size = 1.5) +
  scale_color_manual(values = colors) +
  scale_shape_manual(values = shapes) +
  coord_flip() +
  theme_minimal(base_size = 13) +
  theme(axis.text.y = element_text(size = 20),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 20),
        legend.position = "bottom",
        legend.title = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.major.x = element_line(linewidth = 1.5)
  ) +
  
  labs(
    x = NULL,
    y = NULL
  )+
  theme(legend.position = "none")


dev.off()


#-------------------------------------
# Robustness Checks
#-------------------------------------

#+++++++++++++++++++++++++++++++++++++++++++++++++++
# Ordinal Regression Models  
#+++++++++++++++++++++++++++++++++++++++++++++++++++

#===================================================================================
# Estimations (for Table 8 and Table 9)
#===================================================================================

### Germany
DE_3_cum <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=cumulative(parallel = TRUE), data=DE3)
summary(DE_3_cum)

DE_3_adj <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=acat(parallel = TRUE, reverse = TRUE), data=DE3)
summary(DE_3_adj)

### France
FR_3_cum <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=cumulative(parallel = TRUE), data=FR3)
summary(FR_3_cum)

FR_3_adj <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=acat(parallel = TRUE, reverse = TRUE), data=FR3)
summary(FR_3_adj)

### Italy
IT_3_cum <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=cumulative(parallel = TRUE), data=IT3)
summary(IT_3_cum)

IT_3_adj <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=acat(parallel = TRUE, reverse = TRUE), data=IT3)
summary(IT_3_adj)

### Spain
ES_3_cum <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=cumulative(parallel = TRUE), data=ES3)
summary(ES_3_cum)

ES_3_adj <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=acat(parallel = TRUE, reverse = TRUE), data=ES3)
summary(ES_3_adj)

### Poland
PL_3_cum <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=cumulative(parallel = TRUE), data=PL3)
summary(PL_3_cum)

PL_3_adj <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=acat(parallel = TRUE, reverse = TRUE), data=PL3)
summary(PL_3_adj)

### Romania
RO_3_cum <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=cumulative(parallel = TRUE), data=RO3)
summary(RO_3_cum)

RO_3_adj <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=acat(parallel = TRUE, reverse = TRUE), data=RO3)
summary(RO_3_adj)

### Netherlands
NL_3_cum <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=cumulative(parallel = TRUE), data=NL3)
summary(NL_3_cum)

NL_3_adj <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=acat(parallel = TRUE, reverse = TRUE), data=NL3)
summary(NL_3_adj)

### Sweden
SE_3_cum <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=cumulative(parallel = TRUE), data=SE3)
summary(SE_3_cum)

SE_3_adj <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=acat(parallel = TRUE, reverse = TRUE), data=SE3)
summary(SE_3_adj)

### Austria
AT_3_cum <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=cumulative(parallel = TRUE), data=AT3)
summary(AT_3_cum)

AT_3_adj <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=acat(parallel = TRUE, reverse = TRUE), data=AT3)
summary(AT_3_adj)

### Hungary
HU_3_cum <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=cumulative(parallel = TRUE), data=HU3)
summary(HU_3_cum)

HU_3_adj <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=acat(parallel = TRUE, reverse = TRUE), data=HU3)
summary(HU_3_adj)

### Bulgaria
BG_3_cum <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=cumulative(parallel = TRUE), data=BG3)
summary(BG_3_cum)

BG_3_adj <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=acat(parallel = TRUE, reverse = TRUE), data=BG3)
summary(BG_3_adj)

### Slovakia
SK_3_cum <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=cumulative(parallel = TRUE), data=SK3)
summary(SK_3_cum)

SK_3_adj <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=acat(parallel = TRUE, reverse = TRUE), data=SK3)
summary(SK_3_adj)

### Latvia
LV_3_cum <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=cumulative(parallel = TRUE), data=LV3)
summary(LV_3_cum)

LV_3_adj <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=acat(parallel = TRUE, reverse = TRUE), data=LV3)
summary(LV_3_adj)

### Estonia
EE_3_cum <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=cumulative(parallel = TRUE), data=EE3)
summary(EE_3_cum)

EE_3_adj <- vglm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
                 family=acat(parallel = TRUE, reverse = TRUE), data=EE3)
summary(EE_3_adj)


#---------------------------------------------------------------------------------------
# Table 8: Cumulative Logistic Regression Estimates
#----------------------------------------------------------------------------------------

Table8 <- matrix(NA, 14, ncol = 10)
colnames(Table8) <- c("Coef.", "s.e", "Coef.", "s.e", "Coef.", "s.e", "Coef.", "s.e","Log-Lik", "R^2")
rownames(Table8) <- c("Germany", "France", "Italy", "Spain", "Poland", "Romania", "Netherlands",
                      "Sweden", "Hungary","Austria",  "Bulgaria", "Slovakia", "Latvia", "Estonia")


Table8[1,] <-round(c(coef(DE_3_cum)["pro_imms"],coef(summary(DE_3_cum))[,'Std. Error']["pro_imms"],
                     coef(DE_3_cum)["pro_reds"],coef(summary(DE_3_cum))[,'Std. Error']["pro_reds"],
                     coef(DE_3_cum)["nat_eus"], coef(summary(DE_3_cum))[,'Std. Error']["nat_eus"],
                     coef(DE_3_cum)["nat_regs"],coef(summary(DE_3_cum))[,'Std. Error']["nat_regs"],
                     logLik(DE_3_cum), 1 - (logLik(DE_3_cum) / as.numeric(logLik(update(DE_3_cum, . ~ 1))))), digits = 2)

Table8[2,] <-round(c(coef(FR_3_cum)["pro_imms"],coef(summary(FR_3_cum))[,'Std. Error']["pro_imms"],
                     coef(FR_3_cum)["pro_reds"],coef(summary(FR_3_cum))[,'Std. Error']["pro_reds"] ,
                     coef(FR_3_cum)["nat_eus"], coef(summary(FR_3_cum))[,'Std. Error']["nat_eus"],
                     coef(FR_3_cum)["nat_regs"],coef(summary(FR_3_cum))[,'Std. Error']["nat_regs"] ,
                     logLik(FR_3_cum), 1 - (logLik(FR_3_cum) / as.numeric(logLik(update(FR_3_cum, . ~ 1))))), digits = 2)

Table8[3,] <-round(c(coef(IT_3_cum)["pro_imms"],coef(summary(IT_3_cum))[,'Std. Error']["pro_imms"],
                     coef(IT_3_cum)["pro_reds"],coef(summary(IT_3_cum))[,'Std. Error']["pro_reds"] ,
                     coef(IT_3_cum)["nat_eus"], coef(summary(IT_3_cum))[,'Std. Error']["nat_eus"],
                     coef(IT_3_cum)["nat_regs"],coef(summary(IT_3_cum))[,'Std. Error']["nat_regs"] ,
                     logLik(IT_3_cum), 1 - (logLik(IT_3_cum) / as.numeric(logLik(update(IT_3_cum, . ~ 1))))), digits = 2)

Table8[4,] <-round(c(coef(ES_3_cum)["pro_imms"],coef(summary(ES_3_cum))[,'Std. Error']["pro_imms"],
                     coef(ES_3_cum)["pro_reds"],coef(summary(ES_3_cum))[,'Std. Error']["pro_reds"] ,
                     coef(ES_3_cum)["nat_eus"], coef(summary(ES_3_cum))[,'Std. Error']["nat_eus"],
                     coef(ES_3_cum)["nat_regs"],coef(summary(ES_3_cum))[,'Std. Error']["nat_regs"] ,
                     logLik(ES_3_cum), 1 - (logLik(ES_3_cum) / as.numeric(logLik(update(ES_3_cum, . ~ 1))))), digits = 2)

Table8[5,] <-round(c(coef(PL_3_cum)["pro_imms"],coef(summary(PL_3_cum))[,'Std. Error']["pro_imms"],
                     coef(PL_3_cum)["pro_reds"],coef(summary(PL_3_cum))[,'Std. Error']["pro_reds"] ,
                     coef(PL_3_cum)["nat_eus"], coef(summary(PL_3_cum))[,'Std. Error']["nat_eus"],
                     coef(PL_3_cum)["nat_regs"],coef(summary(PL_3_cum))[,'Std. Error']["nat_regs"] ,
                     logLik(PL_3_cum), 1 - (logLik(PL_3_cum) / as.numeric(logLik(update(PL_3_cum, . ~ 1))))), digits = 2)

Table8[6,] <-round(c(coef(RO_3_cum)["pro_imms"],coef(summary(RO_3_cum))[,'Std. Error']["pro_imms"],
                     coef(RO_3_cum)["pro_reds"],coef(summary(RO_3_cum))[,'Std. Error']["pro_reds"] ,
                     coef(RO_3_cum)["nat_eus"], coef(summary(RO_3_cum))[,'Std. Error']["nat_eus"],
                     coef(RO_3_cum)["nat_regs"],coef(summary(RO_3_cum))[,'Std. Error']["nat_regs"] ,
                     logLik(RO_3_cum), 1 - (logLik(RO_3_cum) / as.numeric(logLik(update(RO_3_cum, . ~ 1))))), digits = 2)

Table8[7,] <-round(c(coef(NL_3_cum)["pro_imms"],coef(summary(NL_3_cum))[,'Std. Error']["pro_imms"],
                     coef(NL_3_cum)["pro_reds"],coef(summary(NL_3_cum))[,'Std. Error']["pro_reds"] ,
                     coef(NL_3_cum)["nat_eus"], coef(summary(NL_3_cum))[,'Std. Error']["nat_eus"],
                     coef(NL_3_cum)["nat_regs"],coef(summary(NL_3_cum))[,'Std. Error']["nat_regs"] ,
                     logLik(NL_3_cum), 1 - (logLik(NL_3_cum) / as.numeric(logLik(update(NL_3_cum, . ~ 1))))), digits = 2)

Table8[8,] <-round(c(coef(SE_3_cum)["pro_imms"],coef(summary(SE_3_cum))[,'Std. Error']["pro_imms"],
                     coef(SE_3_cum)["pro_reds"],coef(summary(SE_3_cum))[,'Std. Error']["pro_reds"] ,
                     coef(SE_3_cum)["nat_eus"], coef(summary(SE_3_cum))[,'Std. Error']["nat_eus"],
                     coef(SE_3_cum)["nat_regs"],coef(summary(SE_3_cum))[,'Std. Error']["nat_regs"] ,
                     logLik(SE_3_cum), 1 - (logLik(SE_3_cum) / as.numeric(logLik(update(SE_3_cum, . ~ 1))))), digits = 2)

Table8[9,] <-round(c(coef(HU_3_cum)["pro_imms"],coef(summary(HU_3_cum))[,'Std. Error']["pro_imms"],
                     coef(HU_3_cum)["pro_reds"],coef(summary(HU_3_cum))[,'Std. Error']["pro_reds"] ,
                     coef(HU_3_cum)["nat_eus"], coef(summary(HU_3_cum))[,'Std. Error']["nat_eus"],
                     coef(HU_3_cum)["nat_regs"],coef(summary(HU_3_cum))[,'Std. Error']["nat_regs"] ,
                     logLik(HU_3_cum), 1 - (logLik(HU_3_cum) / as.numeric(logLik(update(HU_3_cum, . ~ 1))))), digits = 2)

Table8[10,] <-round(c(coef(AT_3_cum)["pro_imms"],coef(summary(AT_3_cum))[,'Std. Error']["pro_imms"],
                      coef(AT_3_cum)["pro_reds"],coef(summary(AT_3_cum))[,'Std. Error']["pro_reds"] ,
                      coef(AT_3_cum)["nat_eus"], coef(summary(AT_3_cum))[,'Std. Error']["nat_eus"],
                      coef(AT_3_cum)["nat_regs"],coef(summary(AT_3_cum))[,'Std. Error']["nat_regs"] ,
                      logLik(AT_3_cum), 1 - (logLik(AT_3_cum) / as.numeric(logLik(update(AT_3_cum, . ~ 1))))), digits = 2)

Table8[11,] <-round(c(coef(BG_3_cum)["pro_imms"],coef(summary(BG_3_cum))[,'Std. Error']["pro_imms"],
                      coef(BG_3_cum)["pro_reds"],coef(summary(BG_3_cum))[,'Std. Error']["pro_reds"] ,
                      coef(BG_3_cum)["nat_eus"], coef(summary(BG_3_cum))[,'Std. Error']["nat_eus"],
                      coef(BG_3_cum)["nat_regs"],coef(summary(BG_3_cum))[,'Std. Error']["nat_regs"] ,
                      logLik(BG_3_cum), 1 - (logLik(BG_3_cum) / as.numeric(logLik(update(BG_3_cum, . ~ 1))))), digits = 2)

Table8[12,] <-round(c(coef(SK_3_cum)["pro_imms"],coef(summary(SK_3_cum))[,'Std. Error']["pro_imms"],
                      coef(SK_3_cum)["pro_reds"],coef(summary(SK_3_cum))[,'Std. Error']["pro_reds"] ,
                      coef(SK_3_cum)["nat_eus"], coef(summary(SK_3_cum))[,'Std. Error']["nat_eus"],
                      coef(SK_3_cum)["nat_regs"],coef(summary(SK_3_cum))[,'Std. Error']["nat_regs"] ,
                      logLik(SK_3_cum), 1 - (logLik(SK_3_cum) / as.numeric(logLik(update(SK_3_cum, . ~ 1))))), digits = 2)

Table8[13,] <-round(c(coef(LV_3_cum)["pro_imms"],coef(summary(LV_3_cum))[,'Std. Error']["pro_imms"],
                      coef(LV_3_cum)["pro_reds"],coef(summary(LV_3_cum))[,'Std. Error']["pro_reds"] ,
                      coef(LV_3_cum)["nat_eus"], coef(summary(LV_3_cum))[,'Std. Error']["nat_eus"],
                      coef(LV_3_cum)["nat_regs"],coef(summary(LV_3_cum))[,'Std. Error']["nat_regs"] ,
                      logLik(LV_3_cum), 1 - (logLik(LV_3_cum) / as.numeric(logLik(update(LV_3_cum, . ~ 1))))), digits = 2)

Table8[14,] <-round(c(coef(EE_3_cum)["pro_imms"],coef(summary(EE_3_cum))[,'Std. Error']["pro_imms"],
                      coef(EE_3_cum)["pro_reds"],coef(summary(EE_3_cum))[,'Std. Error']["pro_reds"] ,
                      coef(EE_3_cum)["nat_eus"], coef(summary(EE_3_cum))[,'Std. Error']["nat_eus"],
                      coef(EE_3_cum)["nat_regs"],coef(summary(EE_3_cum))[,'Std. Error']["nat_regs"] ,
                      logLik(EE_3_cum), 1 - (logLik(EE_3_cum) / as.numeric(logLik(update(EE_3_cum, . ~ 1))))), digits = 2)

Table8


#---------------------------------------------------------------------------------------
# Table 9: Adjacent Categories Regression Estimates
#----------------------------------------------------------------------------------------

Table9 <- matrix(NA, 14, ncol =10)
colnames(Table9) <- c("Coef.", "s.e", "Coef.", "s.e", "Coef.", "s.e", "Coef.", "s.e","Log-Lik", "R^2")
rownames(Table9) <- c("Germany", "France", "Italy", "Spain", "Poland", "Romania", "Netherlands",
                      "Sweden", "Hungary","Austria",  "Bulgaria", "Slovakia", "Latvia", "Estonia")


Table9[1,] <-round(c(coef(DE_3_adj)["pro_imms"],coef(summary(DE_3_adj))[,'Std. Error']["pro_imms"],
                     coef(DE_3_adj)["pro_reds"],coef(summary(DE_3_adj))[,'Std. Error']["pro_reds"] ,
                     coef(DE_3_adj)["nat_eus"], coef(summary(DE_3_adj))[,'Std. Error']["nat_eus"],
                     coef(DE_3_adj)["nat_regs"],coef(summary(DE_3_adj))[,'Std. Error']["nat_regs"] ,
                     logLik(DE_3_adj), 1 - (logLik(DE_3_adj) / as.numeric(logLik(update(DE_3_adj, . ~ 1))))), digits = 2)

Table9[2,] <-round(c(coef(FR_3_adj)["pro_imms"],coef(summary(FR_3_adj))[,'Std. Error']["pro_imms"],
                     coef(FR_3_adj)["pro_reds"],coef(summary(FR_3_adj))[,'Std. Error']["pro_reds"] ,
                     coef(FR_3_adj)["nat_eus"], coef(summary(FR_3_adj))[,'Std. Error']["nat_eus"],
                     coef(FR_3_adj)["nat_regs"],coef(summary(FR_3_adj))[,'Std. Error']["nat_regs"] ,
                     logLik(FR_3_adj), 1 - (logLik(FR_3_adj) / as.numeric(logLik(update(FR_3_adj, . ~ 1))))), digits = 2)

Table9[3,] <-round(c(coef(IT_3_adj)["pro_imms"],coef(summary(IT_3_adj))[,'Std. Error']["pro_imms"],
                     coef(IT_3_adj)["pro_reds"],coef(summary(IT_3_adj))[,'Std. Error']["pro_reds"] ,
                     coef(IT_3_adj)["nat_eus"], coef(summary(IT_3_adj))[,'Std. Error']["nat_eus"],
                     coef(IT_3_adj)["nat_regs"],coef(summary(IT_3_adj))[,'Std. Error']["nat_regs"] ,
                     logLik(IT_3_adj), 1 - (logLik(IT_3_adj) / as.numeric(logLik(update(IT_3_adj, . ~ 1))))), digits = 2)

Table9[4,] <-round(c(coef(ES_3_adj)["pro_imms"],coef(summary(ES_3_adj))[,'Std. Error']["pro_imms"],
                     coef(ES_3_adj)["pro_reds"],coef(summary(ES_3_adj))[,'Std. Error']["pro_reds"] ,
                     coef(ES_3_adj)["nat_eus"], coef(summary(ES_3_adj))[,'Std. Error']["nat_eus"],
                     coef(ES_3_adj)["nat_regs"],coef(summary(ES_3_adj))[,'Std. Error']["nat_regs"] ,
                     logLik(ES_3_adj), 1 - (logLik(ES_3_adj) / as.numeric(logLik(update(ES_3_adj, . ~ 1))))), digits = 2)

Table9[5,] <-round(c(coef(PL_3_adj)["pro_imms"],coef(summary(PL_3_adj))[,'Std. Error']["pro_imms"],
                     coef(PL_3_adj)["pro_reds"],coef(summary(PL_3_adj))[,'Std. Error']["pro_reds"] ,
                     coef(PL_3_adj)["nat_eus"], coef(summary(PL_3_adj))[,'Std. Error']["nat_eus"],
                     coef(PL_3_adj)["nat_regs"],coef(summary(PL_3_adj))[,'Std. Error']["nat_regs"] ,
                     logLik(PL_3_adj), 1 - (logLik(PL_3_adj) / as.numeric(logLik(update(PL_3_adj, . ~ 1))))), digits = 2)

Table9[6,] <-round(c(coef(RO_3_adj)["pro_imms"],coef(summary(RO_3_adj))[,'Std. Error']["pro_imms"],
                     coef(RO_3_adj)["pro_reds"],coef(summary(RO_3_adj))[,'Std. Error']["pro_reds"] ,
                     coef(RO_3_adj)["nat_eus"], coef(summary(RO_3_adj))[,'Std. Error']["nat_eus"],
                     coef(RO_3_adj)["nat_regs"],coef(summary(RO_3_adj))[,'Std. Error']["nat_regs"] ,
                     logLik(RO_3_adj), 1 - (logLik(RO_3_adj) / as.numeric(logLik(update(RO_3_adj, . ~ 1))))), digits = 2)

Table9[7,] <-round(c(coef(NL_3_adj)["pro_imms"],coef(summary(NL_3_adj))[,'Std. Error']["pro_imms"],
                     coef(NL_3_adj)["pro_reds"],coef(summary(NL_3_adj))[,'Std. Error']["pro_reds"] ,
                     coef(NL_3_adj)["nat_eus"], coef(summary(NL_3_adj))[,'Std. Error']["nat_eus"],
                     coef(NL_3_adj)["nat_regs"],coef(summary(NL_3_adj))[,'Std. Error']["nat_regs"] ,
                     logLik(NL_3_adj), 1 - (logLik(NL_3_adj) / as.numeric(logLik(update(NL_3_adj, . ~ 1))))), digits = 2)

Table9[8,] <-round(c(coef(SE_3_adj)["pro_imms"],coef(summary(SE_3_adj))[,'Std. Error']["pro_imms"],
                     coef(SE_3_adj)["pro_reds"],coef(summary(SE_3_adj))[,'Std. Error']["pro_reds"] ,
                     coef(SE_3_adj)["nat_eus"], coef(summary(SE_3_adj))[,'Std. Error']["nat_eus"],
                     coef(SE_3_adj)["nat_regs"],coef(summary(SE_3_adj))[,'Std. Error']["nat_regs"] ,
                     logLik(SE_3_adj), 1 - (logLik(SE_3_adj) / as.numeric(logLik(update(SE_3_adj, . ~ 1))))), digits = 2)

Table9[9,] <-round(c(coef(HU_3_adj)["pro_imms"],coef(summary(HU_3_adj))[,'Std. Error']["pro_imms"],
                     coef(HU_3_adj)["pro_reds"],coef(summary(HU_3_adj))[,'Std. Error']["pro_reds"] ,
                     coef(HU_3_adj)["nat_eus"], coef(summary(HU_3_adj))[,'Std. Error']["nat_eus"],
                     coef(HU_3_adj)["nat_regs"],coef(summary(HU_3_adj))[,'Std. Error']["nat_regs"] ,
                     logLik(HU_3_adj), 1 - (logLik(HU_3_adj) / as.numeric(logLik(update(HU_3_adj, . ~ 1))))), digits = 2)

Table9[10,] <-round(c(coef(AT_3_adj)["pro_imms"],coef(summary(AT_3_adj))[,'Std. Error']["pro_imms"],
                      coef(AT_3_adj)["pro_reds"],coef(summary(AT_3_adj))[,'Std. Error']["pro_reds"] ,
                      coef(AT_3_adj)["nat_eus"], coef(summary(AT_3_adj))[,'Std. Error']["nat_eus"],
                      coef(AT_3_adj)["nat_regs"],coef(summary(AT_3_adj))[,'Std. Error']["nat_regs"] ,
                      logLik(AT_3_adj), 1 - (logLik(AT_3_adj) / as.numeric(logLik(update(AT_3_adj, . ~ 1))))), digits = 2)

Table9[11,] <-round(c(coef(BG_3_adj)["pro_imms"],coef(summary(BG_3_adj))[,'Std. Error']["pro_imms"],
                      coef(BG_3_adj)["pro_reds"],coef(summary(BG_3_adj))[,'Std. Error']["pro_reds"] ,
                      coef(BG_3_adj)["nat_eus"], coef(summary(BG_3_adj))[,'Std. Error']["nat_eus"],
                      coef(BG_3_adj)["nat_regs"],coef(summary(BG_3_adj))[,'Std. Error']["nat_regs"] ,
                      logLik(BG_3_adj), 1 - (logLik(BG_3_adj) / as.numeric(logLik(update(BG_3_adj, . ~ 1))))), digits = 2)

Table9[12,] <-round(c(coef(SK_3_adj)["pro_imms"],coef(summary(SK_3_adj))[,'Std. Error']["pro_imms"],
                      coef(SK_3_adj)["pro_reds"],coef(summary(SK_3_adj))[,'Std. Error']["pro_reds"] ,
                      coef(SK_3_adj)["nat_eus"], coef(summary(SK_3_adj))[,'Std. Error']["nat_eus"],
                      coef(SK_3_adj)["nat_regs"],coef(summary(SK_3_adj))[,'Std. Error']["nat_regs"] ,
                      logLik(SK_3_adj), 1 - (logLik(SK_3_adj) / as.numeric(logLik(update(SK_3_adj, . ~ 1))))), digits = 2)

Table9[13,] <-round(c(coef(LV_3_adj)["pro_imms"],coef(summary(LV_3_adj))[,'Std. Error']["pro_imms"],
                      coef(LV_3_adj)["pro_reds"],coef(summary(LV_3_adj))[,'Std. Error']["pro_reds"] ,
                      coef(LV_3_adj)["nat_eus"], coef(summary(LV_3_adj))[,'Std. Error']["nat_eus"],
                      coef(LV_3_adj)["nat_regs"],coef(summary(LV_3_adj))[,'Std. Error']["nat_regs"] ,
                      logLik(LV_3_adj), 1 - (logLik(LV_3_adj) / as.numeric(logLik(update(LV_3_adj, . ~ 1))))), digits = 2)

Table9[14,] <-round(c(coef(EE_3_adj)["pro_imms"],coef(summary(EE_3_adj))[,'Std. Error']["pro_imms"],
                      coef(EE_3_adj)["pro_reds"],coef(summary(EE_3_adj))[,'Std. Error']["pro_reds"] ,
                      coef(EE_3_adj)["nat_eus"], coef(summary(EE_3_adj))[,'Std. Error']["nat_eus"],
                      coef(EE_3_adj)["nat_regs"],coef(summary(EE_3_adj))[,'Std. Error']["nat_regs"] ,
                      logLik(EE_3_adj), 1 - (logLik(EE_3_adj) / as.numeric(logLik(update(EE_3_adj, . ~ 1))))), digits = 2)

Table9


#---------------------------------------------------------------------------------------
# Estimations with Cohesion Policy (for Table 10)
#----------------------------------------------------------------------------------------

#---------------------------
# Germany (DE)
#---------------------------
DE_y <- glm(eu_support ~ pro_imms + pro_cps + nat_eus + nat_regs  + female + aged + pop + educ + eco + income + eu_ben, data=DEy, family="binomial")

#---------------------------
# France (FR)
#---------------------------
FR_y <- glm(eu_support ~ pro_imms + pro_cps + nat_eus + nat_regs  + female + aged + pop + educ + eco + income + eu_ben, data=FRy, family="binomial")

#---------------------------
# Italy (IT)
#---------------------------
IT_y <- glm(eu_support ~ pro_imms + pro_cps + nat_eus + nat_regs  + female + aged + pop + educ + eco + income + eu_ben, data=ITy, family="binomial")

#---------------------------
# Spain (ES)
#---------------------------
ES_y <- glm(eu_support ~ pro_imms + pro_cps + nat_eus + nat_regs  + female + aged + pop + educ + eco + income + eu_ben, data=ESy, family="binomial")

#---------------------------
# Poland (PL)
#---------------------------
PL_y <- glm(eu_support ~ pro_imms + pro_cps + nat_eus + nat_regs  + female + aged + pop + educ + eco + income + eu_ben, data=PLy, family="binomial")

#---------------------------
# Romania (RO)
#---------------------------
RO_y <- glm(eu_support ~ pro_imms + pro_cps + nat_eus + nat_regs  + female + aged + pop + educ + eco + income + eu_ben, data=ROy, family="binomial")

#---------------------------
# Netherlands (NL)
#---------------------------
NL_y <- glm(eu_support ~ pro_imms + pro_cps + nat_eus + nat_regs  + female + aged + pop + educ + eco + income + eu_ben, data=NLy, family="binomial")

#---------------------------
# Sweden (SE)
#---------------------------
SE_y <- glm(eu_support ~ pro_imms + pro_cps + nat_eus + nat_regs  + female + aged + pop + educ + eco + income + eu_ben, data=SEy, family="binomial")

#---------------------------
# Hungary (HU)
#---------------------------
HU_y <- glm(eu_support ~ pro_imms + pro_cps + nat_eus + nat_regs  + female + aged + pop + educ + eco + income + eu_ben, data=HUy, family="binomial")

#---------------------------
# Austria (AT)
#---------------------------
AT_y <- glm(eu_support ~ pro_imms + pro_cps + nat_eus + nat_regs  + female + aged + pop + educ + eco + income + eu_ben, data=ATy, family="binomial")

#---------------------------
# Bulgaria (BG)
#---------------------------
BG_y <- glm(eu_support ~ pro_imms + pro_cps + nat_eus + nat_regs  + female + aged + pop + educ + eco + income + eu_ben, data=BGy, family="binomial")

#---------------------------
# Slovakia (SK)
#---------------------------
SK_y <- glm(eu_support ~ pro_imms + pro_cps + nat_eus + nat_regs  + female + aged + pop + educ + eco + income + eu_ben, data=SKy, family="binomial")

#---------------------------
# Latvia (LV)
#---------------------------
LV_y <- glm(eu_support ~ pro_imms + pro_cps + nat_eus + nat_regs  + female + aged + pop + educ + eco + income + eu_ben, data=LVy, family="binomial")

#---------------------------
# Estonia (EE)
#---------------------------
EE_y <- glm(eu_support ~ pro_imms + pro_cps + nat_eus + nat_regs  + female + aged + pop + educ + eco + income + eu_ben, data=EEy, family="binomial")


#---------------------------------------------------------------------------------------
# Table 10: Alternative Specification for Redistribution: Agreement with Cohesion Policy
#----------------------------------------------------------------------------------------

Table10 <- matrix(NA, 14, ncol = 10)
colnames(Table10) <- c("beta_1", "", "beta_2", "", "beta_3", "", "beta_4", "", "Log-Lik.", "R^2")
rownames(Table10) <- c("Germany", "France", "Italy", "Spain", "Poland", "Romania", "Netherlands",
                       "Sweden", "Hungary","Austria",  "Bulgaria", "Slovakia", "Latvia", "Estonia")

# DE
Table10[1,] <-round(c(summary(DE_y)$coefficients[2,1], summary(DE_y)$coefficients[2,2],
                      summary(DE_y)$coefficients[3,1], summary(DE_y)$coefficients[3,2],
                      summary(DE_y)$coefficients[4,1], summary(DE_y)$coefficients[4,2],
                      summary(DE_y)$coefficients[5,1], summary(DE_y)$coefficients[5,2],
                      logLik(DE_y), with(summary(DE_y), 1-deviance/null.deviance)),digits=2)
# FR
Table10[2,] <-round(c(summary(FR_y)$coefficients[2,1], summary(FR_y)$coefficients[2,2],
                      summary(FR_y)$coefficients[3,1], summary(FR_y)$coefficients[3,2],
                      summary(FR_y)$coefficients[4,1], summary(FR_y)$coefficients[4,2],
                      summary(FR_y)$coefficients[5,1], summary(FR_y)$coefficients[5,2],
                      logLik(FR_y), with(summary(FR_y), 1-deviance/null.deviance)),digits=2)
# IT
Table10[3,] <-round(c(summary(IT_y)$coefficients[2,1], summary(IT_y)$coefficients[2,2],
                      summary(IT_y)$coefficients[3,1], summary(IT_y)$coefficients[3,2],
                      summary(IT_y)$coefficients[4,1], summary(IT_y)$coefficients[4,2],
                      summary(IT_y)$coefficients[5,1], summary(IT_y)$coefficients[5,2],
                      logLik(IT_y), with(summary(IT_y), 1-deviance/null.deviance)),digits=2)
# ES
Table10[4,] <-round(c(summary(ES_y)$coefficients[2,1], summary(ES_y)$coefficients[2,2],
                      summary(ES_y)$coefficients[3,1], summary(ES_y)$coefficients[3,2],
                      summary(ES_y)$coefficients[4,1], summary(ES_y)$coefficients[4,2],
                      summary(ES_y)$coefficients[5,1], summary(ES_y)$coefficients[5,2],
                      logLik(ES_y), with(summary(ES_y), 1-deviance/null.deviance)),digits=2)
# PL
Table10[5,] <-round(c(summary(PL_y)$coefficients[2,1], summary(PL_y)$coefficients[2,2],
                      summary(PL_y)$coefficients[3,1], summary(PL_y)$coefficients[3,2],
                      summary(PL_y)$coefficients[4,1], summary(PL_y)$coefficients[4,2],
                      summary(PL_y)$coefficients[5,1], summary(PL_y)$coefficients[5,2],
                      logLik(PL_y), with(summary(PL_y), 1-deviance/null.deviance)),digits=2)
# RO
Table10[6,] <-round(c(summary(RO_y)$coefficients[2,1], summary(RO_y)$coefficients[2,2],
                      summary(RO_y)$coefficients[3,1], summary(RO_y)$coefficients[3,2],
                      summary(RO_y)$coefficients[4,1], summary(RO_y)$coefficients[4,2],
                      summary(RO_y)$coefficients[5,1], summary(RO_y)$coefficients[5,2],
                      logLik(RO_y), with(summary(RO_y), 1-deviance/null.deviance)),digits=2)
# NL
Table10[7,] <-round(c(summary(NL_y)$coefficients[2,1], summary(NL_y)$coefficients[2,2],
                      summary(NL_y)$coefficients[3,1], summary(NL_y)$coefficients[3,2],
                      summary(NL_y)$coefficients[4,1], summary(NL_y)$coefficients[4,2],
                      summary(NL_y)$coefficients[5,1], summary(NL_y)$coefficients[5,2],
                      logLik(NL_y), with(summary(NL_y), 1-deviance/null.deviance)),digits=2)
# SE
Table10[8,] <-round(c(summary(SE_y)$coefficients[2,1], summary(SE_y)$coefficients[2,2],
                      summary(SE_y)$coefficients[3,1], summary(SE_y)$coefficients[3,2],
                      summary(SE_y)$coefficients[4,1], summary(SE_y)$coefficients[4,2],
                      summary(SE_y)$coefficients[5,1], summary(SE_y)$coefficients[5,2],
                      logLik(SE_y), with(summary(SE_y), 1-deviance/null.deviance)),digits=2)
# HU
Table10[9,] <-round(c(summary(HU_y)$coefficients[2,1], summary(HU_y)$coefficients[2,2],
                      summary(HU_y)$coefficients[3,1], summary(HU_y)$coefficients[3,2],
                      summary(HU_y)$coefficients[4,1], summary(HU_y)$coefficients[4,2],
                      summary(HU_y)$coefficients[5,1], summary(HU_y)$coefficients[5,2],
                      logLik(HU_y), with(summary(HU_y), 1-deviance/null.deviance)),digits=2)
# AT
Table10[10,] <-round(c(summary(AT_y)$coefficients[2,1], summary(AT_y)$coefficients[2,2],
                       summary(AT_y)$coefficients[3,1], summary(AT_y)$coefficients[3,2],
                       summary(AT_y)$coefficients[4,1], summary(AT_y)$coefficients[4,2],
                       summary(AT_y)$coefficients[5,1], summary(AT_y)$coefficients[5,2],
                       logLik(AT_y), with(summary(AT_y), 1-deviance/null.deviance)),digits=2)
# BG
Table10[11,] <-round(c(summary(BG_y)$coefficients[2,1], summary(BG_y)$coefficients[2,2],
                       summary(BG_y)$coefficients[3,1], summary(BG_y)$coefficients[3,2],
                       summary(BG_y)$coefficients[4,1], summary(BG_y)$coefficients[4,2],
                       summary(BG_y)$coefficients[5,1], summary(BG_y)$coefficients[5,2],
                       logLik(BG_y), with(summary(BG_y), 1-deviance/null.deviance)),digits=2)
# SK
Table10[12,] <-round(c(summary(SK_y)$coefficients[2,1], summary(SK_y)$coefficients[2,2],
                       summary(SK_y)$coefficients[3,1], summary(SK_y)$coefficients[3,2],
                       summary(SK_y)$coefficients[4,1], summary(SK_y)$coefficients[4,2],
                       summary(SK_y)$coefficients[5,1], summary(SK_y)$coefficients[5,2],
                       logLik(SK_y), with(summary(SK_y), 1-deviance/null.deviance)),digits=2)
# LV
Table10[13,] <-round(c(summary(LV_y)$coefficients[2,1], summary(LV_y)$coefficients[2,2],
                       summary(LV_y)$coefficients[3,1], summary(LV_y)$coefficients[3,2],
                       summary(LV_y)$coefficients[4,1], summary(LV_y)$coefficients[4,2],
                       summary(LV_y)$coefficients[5,1], summary(LV_y)$coefficients[5,2],
                       logLik(LV_y), with(summary(LV_y), 1-deviance/null.deviance)),digits=2)
# EE
Table10[14,] <-round(c(summary(EE_y)$coefficients[2,1], summary(EE_y)$coefficients[2,2],
                       summary(EE_y)$coefficients[3,1], summary(EE_y)$coefficients[3,2],
                       summary(EE_y)$coefficients[4,1], summary(EE_y)$coefficients[4,2],
                       summary(EE_y)$coefficients[5,1], summary(EE_y)$coefficients[5,2],
                       logLik(EE_y), with(summary(EE_y), 1-deviance/null.deviance)),digits=2)

Table10


#+++++++++++++++++++++++++++++++++++++++++++++++++++
# ONLINE APPENDIX
#+++++++++++++++++++++++++++++++++++++++++++++++++++

#---------------------------------------------
#  B. Additional Descriptive Statistics
#---------------------------------------------

#---------------------------------------------------------------------------------------
# Table A1: Is EU Membership Good or Bad?
#----------------------------------------------------------------------------------------
TableA1 <- matrix(NA, 14, ncol =5)
colnames(TableA1) <- c("Good", "Bad", "Neither", "Not sure", "Total")
rownames(TableA1) <- c("Germany", "France", "Italy", "Spain", "Poland", "Romania", "Netherlands",
                       "Sweden", "Hungary","Austria",  "Bulgaria", "Slovakia", "Latvia", "Estonia")


TableA1[1,] <-c(summary(DE_A1$eu_support)[1],summary(DE_A1$eu_support)[3],summary(DE_A1$eu_support)[2],
                summary(DE_A1$eu_support)[4],nrow(DE_A1))

TableA1[2,] <-c(summary(FR_A1$eu_support)[1],summary(FR_A1$eu_support)[3],summary(FR_A1$eu_support)[2],
                summary(FR_A1$eu_support)[4],nrow(FR_A1))

TableA1[3,] <-c(summary(IT_A1$eu_support)[1],summary(IT_A1$eu_support)[3],summary(IT_A1$eu_support)[2],
                summary(IT_A1$eu_support)[4],nrow(IT_A1))

TableA1[4,] <-c(summary(ES_A1$eu_support)[1],summary(ES_A1$eu_support)[3],summary(ES_A1$eu_support)[2],
                summary(ES_A1$eu_support)[4],nrow(ES_A1))

TableA1[5,] <-c(summary(PL_A1$eu_support)[1],summary(PL_A1$eu_support)[3],summary(PL_A1$eu_support)[2],
                summary(PL_A1$eu_support)[4],nrow(PL_A1))

TableA1[6,] <-c(summary(RO_A1$eu_support)[1],summary(RO_A1$eu_support)[3],summary(RO_A1$eu_support)[2],
                summary(RO_A1$eu_support)[4],nrow(RO_A1))

TableA1[7,] <-c(summary(NL_A1$eu_support)[1],summary(NL_A1$eu_support)[3],summary(NL_A1$eu_support)[2],
                summary(NL_A1$eu_support)[4],nrow(NL_A1))

TableA1[8,] <-c(summary(SE_A1$eu_support)[1],summary(SE_A1$eu_support)[3],summary(SE_A1$eu_support)[2],
                summary(SE_A1$eu_support)[4],nrow(SE_A1))

TableA1[9,] <-c(summary(HU_A1$eu_support)[1],summary(HU_A1$eu_support)[3],summary(HU_A1$eu_support)[2],
                summary(HU_A1$eu_support)[4],nrow(HU_A1))

TableA1[10,] <-c(summary(AT_A1$eu_support)[1],summary(AT_A1$eu_support)[3],summary(AT_A1$eu_support)[2],
                 summary(AT_A1$eu_support)[4],nrow(AT_A1))

TableA1[11,] <-c(summary(BG_A1$eu_support)[1],summary(BG_A1$eu_support)[3],summary(BG_A1$eu_support)[2],
                 summary(BG_A1$eu_support)[4],nrow(BG_A1))

TableA1[12,] <-c(summary(SK_A1$eu_support)[1],summary(SK_A1$eu_support)[3],summary(SK_A1$eu_support)[2],
                 summary(SK_A1$eu_support)[4],nrow(SK_A1))

TableA1[13,] <-c(summary(LV_A1$eu_support)[1],summary(LV_A1$eu_support)[3],summary(LV_A1$eu_support)[2],
                 summary(LV_A1$eu_support)[4],nrow(LV_A1))

TableA1[14,] <-c(summary(EE_A1$eu_support)[1],summary(EE_A1$eu_support)[3],summary(EE_A1$eu_support)[2],
                 summary(EE_A1$eu_support)[4],nrow(EE_A1))

TableA1


#---------------------------------------------------------------------------------------
# Table A2: Descriptive Statistics on Territorial Identity Differentials
#----------------------------------------------------------------------------------------

# Matrix
TableA2           <- matrix(NA, 14, ncol = 3)
colnames(TableA2) <- c("(n-e)","(n-r)", "corr.")
rownames(TableA2) <- c("Germany", "France", "Italy", "Spain", "Poland", "Romania", "Netherlands", "Sweden", "Hungary", "Austria", "Bulgaria", "Slovakia", "Latvia", "Estonia")

# DE
TableA2[1,]  <-c(round(mean(DE$nat_eu), digits=2),round(mean(DE$nat_reg), digits=2), round(cor(DE$nat_eu,DE$nat_reg), digits=2))
# FR
TableA2[2,]  <-c(round(mean(FR$nat_eu), digits=2),round(mean(FR$nat_reg), digits=2), round(cor(FR$nat_eu,FR$nat_reg), digits=2))
# IT
TableA2[3,]  <-c(round(mean(IT$nat_eu), digits=2),round(mean(IT$nat_reg), digits=2), round(cor(IT$nat_eu,IT$nat_reg), digits=2))
# ES
TableA2[4,]  <-c(round(mean(ES$nat_eu), digits=2),round(mean(ES$nat_reg), digits=2), round(cor(ES$nat_eu,ES$nat_reg), digits=2))
# PL
TableA2[5,]  <-c(round(mean(PL$nat_eu), digits=2),round(mean(PL$nat_reg), digits=2), round(cor(PL$nat_eu,PL$nat_reg), digits=2))
# RO
TableA2[6,]  <-c(round(mean(RO$nat_eu), digits=2),round(mean(RO$nat_reg), digits=2),round(cor(RO$nat_eu,RO$nat_reg), digits=2))
# NL
TableA2[7,]  <-c(round(mean(NL$nat_eu), digits=2),round(mean(NL$nat_reg), digits=2),round(cor(NL$nat_eu,NL$nat_reg), digits=2))
# SE
TableA2[8,]  <-c(round(mean(SE$nat_eu), digits=2),round(mean(SE$nat_reg), digits=2), round(cor(SE$nat_eu,SE$nat_reg), digits=2))
# HU
TableA2[9,]  <-c(round(mean(HU$nat_eu), digits=2),round(mean(HU$nat_reg), digits=2), round(cor(HU$nat_eu,HU$nat_reg), digits=2))
# AT
TableA2[10,]  <-c(round(mean(AT$nat_eu), digits=2),round(mean(AT$nat_reg), digits=2),round(cor(AT$nat_eu,AT$nat_reg), digits=2))
# BG
TableA2[11,]  <-c(round(mean(BG$nat_eu), digits=2),round(mean(BG$nat_reg), digits=2), round(cor(BG$nat_eu,BG$nat_reg), digits=2))
# SK
TableA2[12,]  <-c(round(mean(SK$nat_eu), digits=2),round(mean(SK$nat_reg), digits=2), round(cor(SK$nat_eu,SK$nat_reg), digits=2))
# LV
TableA2[13,]  <-c(round(mean(LV$nat_eu), digits=2),round(mean(LV$nat_reg), digits=2), round(cor(LV$nat_eu,LV$nat_reg), digits=2))
# EE
TableA2[14,]  <-c(round(mean(EE$nat_eu), digits=2),round(mean(EE$nat_reg), digits=2), round(cor(EE$nat_eu,EE$nat_reg), digits=2))

TableA2

#---------------------------------------------------------------------------------------
# Table A3: Association Between Level of Income and Preference for Redistribution
#----------------------------------------------------------------------------------------

TableA3 <- matrix(NA, 14, ncol = 6)
colnames(TableA3) <- c("Low","sd","Medium","sd","High","sd")
rownames(TableA3) <- c("Germany", "France", "Italy", "Spain", "Poland", "Romania", "Netherlands", "Sweden", "Hungary", "Austria", "Bulgaria", "Slovakia", "Latvia", "Estonia")

# DE
TableA3[1,] <-round(c(mean(DE$pro_red[DE$income=="Low"]),sd(DE$pro_red[DE$income=="Low"]),
                      mean(DE$pro_red[DE$income=="Medium"]),sd(DE$pro_red[DE$income=="Medium"]),
                      mean(DE$pro_red[DE$income=="High"]),sd(DE$pro_red[DE$income=="High"])), digits=2)

# FR
TableA3[2,]  <-round(c(mean(FR$pro_red[FR$income=="Low"]),sd(FR$pro_red[FR$income=="Low"]),
                       mean(FR$pro_red[FR$income=="Medium"]),sd(FR$pro_red[FR$income=="Medium"]),
                       mean(FR$pro_red[FR$income=="High"]),sd(FR$pro_red[FR$income=="High"])), digits=2)

# IT
TableA3[3,]  <-round(c(mean(IT$pro_red[IT$income=="Low"]),sd(IT$pro_red[IT$income=="Low"]),
                       mean(IT$pro_red[IT$income=="Medium"]),sd(IT$pro_red[IT$income=="Medium"]),
                       mean(IT$pro_red[IT$income=="High"]),sd(IT$pro_red[IT$income=="High"])), digits=2)

# ES
TableA3[4,] <-round(c(mean(ES$pro_red[ES$income=="Low"]),sd(ES$pro_red[ES$income=="Low"]),
                      mean(ES$pro_red[ES$income=="Medium"]),sd(ES$pro_red[ES$income=="Medium"]),
                      mean(ES$pro_red[ES$income=="High"]),sd(ES$pro_red[ES$income=="High"])), digits=2)

# PL
TableA3[5,] <-round(c(mean(PL$pro_red[PL$income=="Low"]),sd(PL$pro_red[PL$income=="Low"]),
                      mean(PL$pro_red[PL$income=="Medium"]),sd(PL$pro_red[PL$income=="Medium"]),
                      mean(PL$pro_red[PL$income=="High"]),sd(PL$pro_red[PL$income=="High"])), digits=2)

# RO
TableA3[6,] <-round(c(mean(RO$pro_red[RO$income=="Low"]),sd(RO$pro_red[RO$income=="Low"]),
                      mean(RO$pro_red[RO$income=="Medium"]),sd(RO$pro_red[RO$income=="Medium"]),
                      mean(RO$pro_red[RO$income=="High"]),sd(RO$pro_red[RO$income=="High"])), digits=2)

# NL
TableA3[7,]<-round(c(mean(NL$pro_red[NL$income=="Low"]),sd(NL$pro_red[NL$income=="Low"]),
                     mean(NL$pro_red[NL$income=="Medium"]),sd(NL$pro_red[NL$income=="Medium"]),
                     mean(NL$pro_red[NL$income=="High"]),sd(NL$pro_red[NL$income=="High"])), digits=2)

# SE
TableA3[8,] <-round(c(mean(SE$pro_red[SE$income=="Low"]),sd(SE$pro_red[SE$income=="Low"]),
                      mean(SE$pro_red[SE$income=="Medium"]),sd(SE$pro_red[SE$income=="Medium"]),
                      mean(SE$pro_red[SE$income=="High"]),sd(SE$pro_red[SE$income=="High"])), digits=2)

# HU
TableA3[9,] <-round(c(mean(HU$pro_red[HU$income=="Low"]),sd(HU$pro_red[HU$income=="Low"]),
                      mean(HU$pro_red[HU$income=="Medium"]),sd(HU$pro_red[HU$income=="Medium"]),
                      mean(HU$pro_red[HU$income=="High"]),sd(HU$pro_red[HU$income=="High"])), digits=2)

# AT
TableA3[10,] <-round(c(mean(AT$pro_red[AT$income=="Low"]),sd(AT$pro_red[AT$income=="Low"]),
                       mean(AT$pro_red[AT$income=="Medium"]),sd(AT$pro_red[AT$income=="Medium"]),
                       mean(AT$pro_red[AT$income=="High"]),sd(AT$pro_red[AT$income=="High"])), digits=2)

# BG
TableA3[11,] <-round(c(mean(BG$pro_red[BG$income=="Low"]),sd(BG$pro_red[BG$income=="Low"]),
                       mean(BG$pro_red[BG$income=="Medium"]),sd(BG$pro_red[BG$income=="Medium"]),
                       mean(BG$pro_red[BG$income=="High"]),sd(BG$pro_red[BG$income=="High"])), digits=2)

# SK
TableA3[12,] <-round(c(mean(SK$pro_red[SK$income=="Low"]),sd(SK$pro_red[SK$income=="Low"]),
                       mean(SK$pro_red[SK$income=="Medium"]),sd(SK$pro_red[SK$income=="Medium"]),
                       mean(SK$pro_red[SK$income=="High"]),sd(SK$pro_red[SK$income=="High"])), digits=2)

# LV
TableA3[13,] <-round(c(mean(LV$pro_red[LV$income=="Low"]),sd(LV$pro_red[LV$income=="Low"]),
                       mean(LV$pro_red[LV$income=="Medium"]),sd(LV$pro_red[LV$income=="Medium"]),
                       mean(LV$pro_red[LV$income=="High"]),sd(LV$pro_red[LV$income=="High"])), digits=2)

# EE
TableA3[14,] <-round(c(mean(EE$pro_red[EE$income=="Low"]),sd(EE$pro_red[EE$income=="Low"]),
                       mean(EE$pro_red[EE$income=="Medium"]),sd(EE$pro_red[EE$income=="Medium"]),
                       mean(EE$pro_red[EE$income=="High"]),sd(EE$pro_red[EE$income=="High"])), digits=2)

TableA3


#------------------------------
# C. Estimation Tables
#------------------------------

#------------------------------
# C.1 Model Comparisons
#------------------------------

#---------------------------------------------------------------------------------------
# Table A4: Explaining Support for EU Membership: Model Comparisons based on BIC
#----------------------------------------------------------------------------------------

# BIC
TableA4 <- matrix(NA, 14, ncol = 4)
colnames(TableA4) <- c("obs","M1: Only controls", "M2: M1 + x_i + y_i", "M3: M2 + (n-e) + (n-r)")
rownames(TableA4) <- c("Germany", "France", "Italy", "Spain", "Poland", "Romania", "Netherlands",
                       "Sweden", "Hungary","Austria",  "Bulgaria", "Slovakia", "Latvia", "Estonia")
TableA4[1,]   <- c(nrow(DE),round(c(BIC(DE_1), BIC(DE_2), BIC(DE_3)), digits=2))
TableA4[2,]   <- c(nrow(FR),round(c(BIC(FR_1), BIC(FR_2), BIC(FR_3)), digits=2))
TableA4[3,]   <- c(nrow(IT),round(c(BIC(IT_1), BIC(IT_2), BIC(IT_3)), digits=2))
TableA4[4,]   <- c(nrow(ES),round(c(BIC(ES_1), BIC(ES_2), BIC(ES_3)), digits=2))
TableA4[5,]   <- c(nrow(PL),round(c(BIC(PL_1), BIC(PL_2), BIC(PL_3)), digits=2))
TableA4[6,]   <- c(nrow(RO),round(c(BIC(RO_1), BIC(RO_2), BIC(RO_3)), digits=2))
TableA4[7,]   <- c(nrow(NL),round(c(BIC(NL_1), BIC(NL_2), BIC(NL_3)), digits=2))
TableA4[8,]   <- c(nrow(SE),round(c(BIC(SE_1), BIC(SE_2), BIC(SE_3)), digits=2))
TableA4[9,]   <- c(nrow(HU),round(c(BIC(HU_1), BIC(HU_2), BIC(HU_3)), digits=2))
TableA4[10,]  <- c(nrow(AT),round(c(BIC(AT_1), BIC(AT_2), BIC(AT_3)), digits=2))
TableA4[11,]  <- c(nrow(BG),round(c(BIC(BG_1), BIC(BG_2), BIC(BG_3)), digits=2))
TableA4[12,]  <- c(nrow(SK),round(c(BIC(SK_1), BIC(SK_2), BIC(SK_3)), digits=2))
TableA4[13,]  <- c(nrow(LV),round(c(BIC(LV_1), BIC(LV_2), BIC(LV_3)), digits=2))
TableA4[14,]  <- c(nrow(EE),round(c(BIC(EE_1), BIC(EE_2), BIC(EE_3)), digits=2))
TableA4


##++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# C.2 Detailed Logistic Regressions by Country
##++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# M3

#----------------
# Germany (DE)
#----------------
TableA5 <- matrix(NA, 15, ncol = 3)
colnames(TableA5) <- c("coef.", "s.e.", "p-value" )
rownames(TableA5) <- c( "Intercept", "Immigration", "Redistribution", "National-European", "National-Regional",
                        "Female", "Age", "Population Size: small city", "Population Size: large city",
                        "Education: high school", "Education: university","Economic Situation: satisfied", 
                        "Income: medium", "Income: high", "Benefit EU Funds: yes")
TableA5[,1]   <- round(summary(DE_3)$coefficients[,1], digits=2)
TableA5[,2]   <- round(summary(DE_3)$coefficients[,2], digits=2)
TableA5[,3]   <- round(summary(DE_3)$coefficients[,4], digits=2)
TableA5

#-----------------
# France (FR)
#----------------
TableA6 <- matrix(NA, 15, ncol = 3)
colnames(TableA6) <- c("coef.", "s.e.", "p-value" )
rownames(TableA6) <- c( "Intercept", "Immigration", "Redistribution", "National-European", "National-Regional",
                        "Female", "Age", "Population Size: small city", "Population Size: large city",
                        "Education: high school", "Education: university","Economic Situation: satisfied",
                        "Income: medium", "Income: high", "Benefit EU Funds: yes")
TableA6[,1]   <- round(summary(FR_3)$coefficients[,1], digits=2)
TableA6[,2]   <- round(summary(FR_3)$coefficients[,2], digits=2)
TableA6[,3]   <- round(summary(FR_3)$coefficients[,4], digits=2)
TableA6

#----------------
# Italy (IT)
#----------------
TableA7 <- matrix(NA, 15, ncol = 3)
colnames(TableA7) <- c("coef.", "s.e.", "p-value" )
rownames(TableA7) <- c( "Intercept", "Immigration", "Redistribution", "National-European", "National-Regional",
                        "Female", "Age", "Population Size: small city", "Population Size: large city",
                        "Education: high school", "Education: university","Economic Situation: satisfied",
                        "Income: medium", "Income: high", "Benefit EU Funds: yes")
TableA7[,1]   <- round(summary(IT_3)$coefficients[,1], digits=2)
TableA7[,2]   <- round(summary(IT_3)$coefficients[,2], digits=2)
TableA7[,3]   <- round(summary(IT_3)$coefficients[,4], digits=2)
TableA7

#----------------
# Spain (ES)
#----------------
TableA8 <- matrix(NA, 15, ncol = 3)
colnames(TableA8) <- c("coef.", "s.e.", "p-value" )
rownames(TableA8) <- c( "Intercept", "Immigration", "Redistribution", "National-European", "National-Regional",
                        "Female", "Age", "Population Size: small city", "Population Size: large city",
                        "Education: high school", "Education: university","Economic Situation: satisfied",
                        "Income: medium", "Income: high", "Benefit EU Funds: yes")
TableA8[,1]   <- round(summary(ES_3)$coefficients[,1], digits=2)
TableA8[,2]   <- round(summary(ES_3)$coefficients[,2], digits=2)
TableA8[,3]   <- round(summary(ES_3)$coefficients[,4], digits=2)
TableA8

#----------------
# Poland (PL)
#----------------
TableA9 <- matrix(NA, 15, ncol = 3)
colnames(TableA9) <- c("coef.", "s.e.", "p-value" )
rownames(TableA9) <- c( "Intercept", "Immigration", "Redistribution", "National-European", "National-Regional",
                        "Female", "Age", "Population Size: small city", "Population Size: large city",
                        "Education: high school", "Education: university","Economic Situation: satisfied",
                        "Income: medium", "Income: high", "Benefit EU Funds: yes")
TableA9[,1]   <- round(summary(PL_3)$coefficients[,1], digits=2)
TableA9[,2]   <- round(summary(PL_3)$coefficients[,2], digits=2)
TableA9[,3]   <- round(summary(PL_3)$coefficients[,4], digits=2)
TableA9

#----------------
# Romania (RO)
#----------------
TableA10 <- matrix(NA, 15, ncol = 3)
colnames(TableA10) <- c("coef.", "s.e.", "p-value" )
rownames(TableA10) <- c( "Intercept", "Immigration", "Redistribution", "National-European", "National-Regional",
                         "Female", "Age", "Population Size: small city", "Population Size: large city",
                         "Education: high school", "Education: university","Economic Situation: satisfied",
                         "Income: medium", "Income: high", "Benefit EU Funds: yes")
TableA10[,1]   <- round(summary(RO_3)$coefficients[,1], digits=2)
TableA10[,2]   <- round(summary(RO_3)$coefficients[,2], digits=2)
TableA10[,3]   <- round(summary(RO_3)$coefficients[,4], digits=2)
TableA10

#----------------
# Netherlands (NL)
#----------------
TableA11 <- matrix(NA, 15, ncol = 3)
colnames(TableA11) <- c("coef.", "s.e.", "p-value" )
rownames(TableA11) <- c( "Intercept", "Immigration", "Redistribution", "National-European", "National-Regional",
                         "Female", "Age", "Population Size: small city", "Population Size: large city",
                         "Education: high school", "Education: university","Economic Situation: satisfied",
                         "Income: medium", "Income: high", "Benefit EU Funds: yes")
TableA11[,1]   <- round(summary(NL_3)$coefficients[,1], digits=2)
TableA11[,2]   <- round(summary(NL_3)$coefficients[,2], digits=2)
TableA11[,3]   <- round(summary(NL_3)$coefficients[,4], digits=2)
TableA11

#----------------
# Sweden (SE)
#----------------
TableA12 <- matrix(NA, 15, ncol = 3)
colnames(TableA12) <- c("coef.", "s.e.", "p-value" )
rownames(TableA12) <- c( "Intercept", "Immigration", "Redistribution", "National-European", "National-Regional",
                         "Female", "Age", "Population Size: small city", "Population Size: large city",
                         "Education: high school", "Education: university","Economic Situation: satisfied",
                         "Income: medium", "Income: high", "Benefit EU Funds: yes")
TableA12[,1]   <- round(summary(SE_3)$coefficients[,1], digits=2)
TableA12[,2]   <- round(summary(SE_3)$coefficients[,2], digits=2)
TableA12[,3]   <- round(summary(SE_3)$coefficients[,4], digits=2)
TableA12

#----------------
# Hungary (HU)
#----------------
TableA13 <- matrix(NA, 15, ncol = 3)
colnames(TableA13) <- c("coef.", "s.e.", "p-value" )
rownames(TableA13) <- c( "Intercept", "Immigration", "Redistribution", "National-European", "National-Regional",
                         "Female", "Age", "Population Size: small city", "Population Size: large city",
                         "Education: high school", "Education: university","Economic Situation: satisfied",
                         "Income: medium", "Income: high", "Benefit EU Funds: yes")
TableA13[,1]   <- round(summary(HU_3)$coefficients[,1], digits=2)
TableA13[,2]   <- round(summary(HU_3)$coefficients[,2], digits=2)
TableA13[,3]   <- round(summary(HU_3)$coefficients[,4], digits=2)
TableA13


#----------------
# Austria (AT)
#----------------
TableA14 <- matrix(NA, 15, ncol = 3)
colnames(TableA14) <- c("coef.", "s.e.", "p-value" )
rownames(TableA14) <- c( "Intercept", "Immigration", "Redistribution", "National-European", "National-Regional",
                         "Female", "Age", "Population Size: small city", "Population Size: large city",
                         "Education: high school", "Education: university","Economic Situation: satisfied",
                         "Income: medium", "Income: high", "Benefit EU Funds: yes")
TableA14[,1]   <- round(summary(AT_3)$coefficients[,1], digits=2)
TableA14[,2]   <- round(summary(AT_3)$coefficients[,2], digits=2)
TableA14[,3]   <- round(summary(AT_3)$coefficients[,4], digits=2)
TableA14

#----------------
# Bulgaria (BG)
#----------------
TableA15 <- matrix(NA, 15, ncol = 3)
colnames(TableA15) <- c("coef.", "s.e.", "p-value" )
rownames(TableA15) <- c( "Intercept", "Immigration", "Redistribution", "National-European", "National-Regional",
                         "Female", "Age", "Population Size: small city", "Population Size: large city",
                         "Education: high school", "Education: university","Economic Situation: satisfied",
                         "Income: medium", "Income: high", "Benefit EU Funds: yes")
TableA15[,1]   <- round(summary(BG_3)$coefficients[,1], digits=2)
TableA15[,2]   <- round(summary(BG_3)$coefficients[,2], digits=2)
TableA15[,3]   <- round(summary(BG_3)$coefficients[,4], digits=2)
TableA15

#----------------
# Slovakia (SK)
#----------------
TableA16 <- matrix(NA, 15, ncol = 3)
colnames(TableA16) <- c("coef.", "s.e.", "p-value" )
rownames(TableA16) <- c( "Intercept", "Immigration", "Redistribution", "National-European", "National-Regional",
                         "Female", "Age", "Population Size: small city", "Population Size: large city",
                         "Education: high school", "Education: university","Economic Situation: satisfied",
                         "Income: medium", "Income: high", "Benefit EU Funds: yes")
TableA16[,1]   <- round(summary(SK_3)$coefficients[,1], digits=2)
TableA16[,2]   <- round(summary(SK_3)$coefficients[,2], digits=2)
TableA16[,3]   <- round(summary(SK_3)$coefficients[,4], digits=2)
TableA16

#---------------
# Latvia (LV)
#----------------
TableA17 <- matrix(NA, 15, ncol = 3)
colnames(TableA17) <- c("coef.", "s.e.", "p-value" )
rownames(TableA17) <- c( "Intercept", "Immigration", "Redistribution", "National-European", "National-Regional",
                         "Female", "Age", "Population Size: small city", "Population Size: large city",
                         "Education: high school", "Education: university","Economic Situation: satisfied",
                         "Income: medium", "Income: high", "Benefit EU Funds: yes")
TableA17[,1]   <- round(summary(LV_3)$coefficients[,1], digits=2)
TableA17[,2]   <- round(summary(LV_3)$coefficients[,2], digits=2)
TableA17[,3]   <- round(summary(LV_3)$coefficients[,4], digits=2)
TableA17

#----------------
# Estonia (EE)
#----------------
TableA18 <- matrix(NA, 15, ncol = 3)
colnames(TableA18) <- c("coef.", "s.e.", "p-value" )
rownames(TableA18) <- c( "Intercept", "Immigration", "Redistribution", "National-European", "National-Regional",
                         "Female", "Age", "Population Size: small city", "Population Size: large city",
                         "Education: high school", "Education: university","Economic Situation: satisfied",
                         "Income: medium", "Income: high", "Benefit EU Funds: yes")
TableA18[,1]   <- round(summary(EE_3)$coefficients[,1], digits=2)
TableA18[,2]   <- round(summary(EE_3)$coefficients[,2], digits=2)
TableA18[,3]   <- round(summary(EE_3)$coefficients[,4], digits=2)
TableA18


##++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# C.3 Estimates with All Identities Included
##++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#---------------------------------------------------------------------------------------
# Estimations with Three Identities (for Table A19)
#----------------------------------------------------------------------------------------

#---------------------------
# Germany (DE)
#---------------------------
DE_3id <- glm(eu_support ~ pro_imm + pro_red + nat_id + reg_id + eur_id + female + aged + pop + educ + eco + income + eu_ben, data=DE, family="binomial")

#---------------------------
# France (FR)
#---------------------------
FR_3id <- glm(eu_support ~ pro_imm + pro_red + nat_id + reg_id + eur_id + female + aged + pop + educ + eco + income + eu_ben, data=FR, family="binomial")

#---------------------------
# Italy (IT)
#---------------------------
IT_3id <- glm(eu_support ~ pro_imm + pro_red + nat_id + reg_id + eur_id + female + aged + pop + educ + eco + income + eu_ben, data=IT, family="binomial")

#---------------------------
# Spain (ES)
#---------------------------
ES_3id <- glm(eu_support ~ pro_imm + pro_red + nat_id + reg_id + eur_id + female + aged + pop + educ + eco + income + eu_ben, data=ES, family="binomial")

#---------------------------
# Poland (PL)
#---------------------------
PL_3id <- glm(eu_support ~ pro_imm + pro_red + nat_id + reg_id + eur_id + female + aged + pop + educ + eco + income + eu_ben, data=PL, family="binomial")

#---------------------------
# Romania (RO)
#---------------------------
RO_3id <- glm(eu_support ~ pro_imm + pro_red + nat_id + reg_id + eur_id + female + aged + pop + educ + eco + income + eu_ben, data=RO, family="binomial")

#---------------------------
# Netherlands (NL)
#---------------------------
NL_3id <- glm(eu_support ~ pro_imm + pro_red + nat_id + reg_id + eur_id + female + aged + pop + educ + eco + income + eu_ben, data=NL, family="binomial")

#---------------------------
# Sweden (SE)
#---------------------------
SE_3id <- glm(eu_support ~ pro_imm + pro_red + nat_id + reg_id + eur_id + female + aged + pop + educ + eco + income + eu_ben, data=SE, family="binomial")

#---------------------------
# Hungary (HU)
#---------------------------
HU_3id <- glm(eu_support ~ pro_imm + pro_red + nat_id + reg_id + eur_id + female + aged + pop + educ + eco + income + eu_ben, data=HU, family="binomial")

#---------------------------
# Austria (AT)
#---------------------------
AT_3id <- glm(eu_support ~ pro_imm + pro_red + nat_id + reg_id + eur_id + female + aged + pop + educ + eco + income + eu_ben, data=AT, family="binomial")

#---------------------------
# Bulgaria (BG)
#---------------------------
BG_3id <- glm(eu_support ~ pro_imm + pro_red + nat_id + reg_id + eur_id + female + aged + pop + educ + eco + income + eu_ben, data=BG, family="binomial")

#---------------------------
# Slovakia (SK)
#---------------------------
SK_3id <- glm(eu_support ~ pro_imm + pro_red + nat_id + reg_id + eur_id + female + aged + pop + educ + eco + income + eu_ben, data=SK, family="binomial")

#---------------------------
# Latvia (LV)
#---------------------------
LV_3id <- glm(eu_support ~ pro_imm + pro_red + nat_id + reg_id + eur_id + female + aged + pop + educ + eco + income + eu_ben, data=LV, family="binomial")

#---------------------------
# Estonia (EE)
#---------------------------
EE_3id <- glm(eu_support ~ pro_imm + pro_red + nat_id + reg_id + eur_id + female + aged + pop + educ + eco + income + eu_ben, data=EE, family="binomial")


#---------------------------------------------------------------------------------------
# Table A19: Binary Logistic Regression Models: Three Identities
#----------------------------------------------------------------------------------------

TableA19 <- matrix(NA, 14, ncol = 12)
colnames(TableA19) <- c("beta_1", "", "beta_2", "", "beta_3", "", "beta_4", "", "beta_5","", "Log-Lik.", "R^2")
rownames(TableA19) <- c("Germany", "France", "Italy", "Spain", "Poland", "Romania", "Netherlands",
                        "Sweden", "Hungary","Austria",  "Bulgaria", "Slovakia", "Latvia", "Estonia")

# DE
TableA19[1,] <-round(c(summary(DE_3id)$coefficients[2,1], summary(DE_3id)$coefficients[2,2],
                       summary(DE_3id)$coefficients[3,1], summary(DE_3id)$coefficients[3,2],
                       summary(DE_3id)$coefficients[4,1], summary(DE_3id)$coefficients[4,2],
                       summary(DE_3id)$coefficients[5,1], summary(DE_3id)$coefficients[5,2],
                       summary(DE_3id)$coefficients[6,1], summary(DE_3id)$coefficients[6,2],
                       logLik(DE_3id), with(summary(DE_3id), 1-deviance/null.deviance)),digits=2)
# FR
TableA19[2,] <-round(c(summary(FR_3id)$coefficients[2,1], summary(FR_3id)$coefficients[2,2],
                       summary(FR_3id)$coefficients[3,1], summary(FR_3id)$coefficients[3,2],
                       summary(FR_3id)$coefficients[4,1], summary(FR_3id)$coefficients[4,2],
                       summary(FR_3id)$coefficients[5,1], summary(FR_3id)$coefficients[5,2],
                       summary(FR_3id)$coefficients[6,1], summary(FR_3id)$coefficients[6,2],
                       logLik(FR_3id), with(summary(FR_3id), 1-deviance/null.deviance)),digits=2)
# IT
TableA19[3,] <-round(c(summary(IT_3id)$coefficients[2,1], summary(IT_3id)$coefficients[2,2],
                       summary(IT_3id)$coefficients[3,1], summary(IT_3id)$coefficients[3,2],
                       summary(IT_3id)$coefficients[4,1], summary(IT_3id)$coefficients[4,2],
                       summary(IT_3id)$coefficients[5,1], summary(IT_3id)$coefficients[5,2],
                       summary(IT_3id)$coefficients[6,1], summary(IT_3id)$coefficients[6,2],
                       logLik(IT_3id), with(summary(IT_3id), 1-deviance/null.deviance)),digits=2)
# ES
TableA19[4,] <-round(c(summary(ES_3id)$coefficients[2,1], summary(ES_3id)$coefficients[2,2],
                       summary(ES_3id)$coefficients[3,1], summary(ES_3id)$coefficients[3,2],
                       summary(ES_3id)$coefficients[4,1], summary(ES_3id)$coefficients[4,2],
                       summary(ES_3id)$coefficients[5,1], summary(ES_3id)$coefficients[5,2],
                       summary(ES_3id)$coefficients[6,1], summary(ES_3id)$coefficients[6,2],
                       logLik(ES_3id), with(summary(ES_3id), 1-deviance/null.deviance)),digits=2)
# PL
TableA19[5,] <-round(c(summary(PL_3id)$coefficients[2,1], summary(PL_3id)$coefficients[2,2],
                       summary(PL_3id)$coefficients[3,1], summary(PL_3id)$coefficients[3,2],
                       summary(PL_3id)$coefficients[4,1], summary(PL_3id)$coefficients[4,2],
                       summary(PL_3id)$coefficients[5,1], summary(PL_3id)$coefficients[5,2],
                       summary(PL_3id)$coefficients[6,1], summary(PL_3id)$coefficients[6,2],
                       logLik(PL_3id), with(summary(PL_3id), 1-deviance/null.deviance)),digits=2)
# RO
TableA19[6,] <-round(c(summary(RO_3id)$coefficients[2,1], summary(RO_3id)$coefficients[2,2],
                       summary(RO_3id)$coefficients[3,1], summary(RO_3id)$coefficients[3,2],
                       summary(RO_3id)$coefficients[4,1], summary(RO_3id)$coefficients[4,2],
                       summary(RO_3id)$coefficients[5,1], summary(RO_3id)$coefficients[5,2],
                       summary(RO_3id)$coefficients[6,1], summary(RO_3id)$coefficients[6,2],
                       logLik(RO_3id), with(summary(RO_3id), 1-deviance/null.deviance)),digits=2)
# NL
TableA19[7,] <-round(c(summary(NL_3id)$coefficients[2,1], summary(NL_3id)$coefficients[2,2],
                       summary(NL_3id)$coefficients[3,1], summary(NL_3id)$coefficients[3,2],
                       summary(NL_3id)$coefficients[4,1], summary(NL_3id)$coefficients[4,2],
                       summary(NL_3id)$coefficients[5,1], summary(NL_3id)$coefficients[5,2],
                       summary(NL_3id)$coefficients[6,1], summary(NL_3id)$coefficients[6,2],
                       logLik(NL_3id), with(summary(NL_3id), 1-deviance/null.deviance)),digits=2)
# SE
TableA19[8,] <-round(c(summary(SE_3id)$coefficients[2,1], summary(SE_3id)$coefficients[2,2],
                       summary(SE_3id)$coefficients[3,1], summary(SE_3id)$coefficients[3,2],
                       summary(SE_3id)$coefficients[4,1], summary(SE_3id)$coefficients[4,2],
                       summary(SE_3id)$coefficients[5,1], summary(SE_3id)$coefficients[5,2],
                       summary(SE_3id)$coefficients[6,1], summary(SE_3id)$coefficients[6,2],
                       logLik(SE_3id), with(summary(SE_3id), 1-deviance/null.deviance)),digits=2)
# HU
TableA19[9,] <-round(c(summary(HU_3id)$coefficients[2,1], summary(HU_3id)$coefficients[2,2],
                       summary(HU_3id)$coefficients[3,1], summary(HU_3id)$coefficients[3,2],
                       summary(HU_3id)$coefficients[4,1], summary(HU_3id)$coefficients[4,2],
                       summary(HU_3id)$coefficients[5,1], summary(HU_3id)$coefficients[5,2],
                       summary(HU_3id)$coefficients[6,1], summary(HU_3id)$coefficients[6,2],
                       logLik(HU_3id), with(summary(HU_3id), 1-deviance/null.deviance)),digits=2)
# AT
TableA19[10,] <-round(c(summary(AT_3id)$coefficients[2,1], summary(AT_3id)$coefficients[2,2],
                        summary(AT_3id)$coefficients[3,1], summary(AT_3id)$coefficients[3,2],
                        summary(AT_3id)$coefficients[4,1], summary(AT_3id)$coefficients[4,2],
                        summary(AT_3id)$coefficients[5,1], summary(AT_3id)$coefficients[5,2],
                        summary(AT_3id)$coefficients[6,1], summary(AT_3id)$coefficients[6,2],
                        logLik(AT_3id), with(summary(AT_3id), 1-deviance/null.deviance)),digits=2)
# BG
TableA19[11,] <-round(c(summary(BG_3id)$coefficients[2,1], summary(BG_3id)$coefficients[2,2],
                        summary(BG_3id)$coefficients[3,1], summary(BG_3id)$coefficients[3,2],
                        summary(BG_3id)$coefficients[4,1], summary(BG_3id)$coefficients[4,2],
                        summary(BG_3id)$coefficients[5,1], summary(BG_3id)$coefficients[5,2],
                        summary(BG_3id)$coefficients[6,1], summary(BG_3id)$coefficients[6,2],
                        logLik(BG_3id), with(summary(BG_3id), 1-deviance/null.deviance)),digits=2)
# SK
TableA19[12,] <-round(c(summary(SK_3id)$coefficients[2,1], summary(SK_3id)$coefficients[2,2],
                        summary(SK_3id)$coefficients[3,1], summary(SK_3id)$coefficients[3,2],
                        summary(SK_3id)$coefficients[4,1], summary(SK_3id)$coefficients[4,2],
                        summary(SK_3id)$coefficients[5,1], summary(SK_3id)$coefficients[5,2],
                        summary(SK_3id)$coefficients[6,1], summary(SK_3id)$coefficients[6,2],
                        logLik(SK_3id), with(summary(SK_3id), 1-deviance/null.deviance)),digits=2)
# LV
TableA19[13,] <-round(c(summary(LV_3id)$coefficients[2,1], summary(LV_3id)$coefficients[2,2],
                        summary(LV_3id)$coefficients[3,1], summary(LV_3id)$coefficients[3,2],
                        summary(LV_3id)$coefficients[4,1], summary(LV_3id)$coefficients[4,2],
                        summary(LV_3id)$coefficients[5,1], summary(LV_3id)$coefficients[5,2],
                        summary(LV_3id)$coefficients[6,1], summary(LV_3id)$coefficients[6,2],
                        logLik(LV_3id), with(summary(LV_3id), 1-deviance/null.deviance)),digits=2)
# EE
TableA19[14,] <-round(c(summary(EE_3id)$coefficients[2,1], summary(EE_3id)$coefficients[2,2],
                        summary(EE_3id)$coefficients[3,1], summary(EE_3id)$coefficients[3,2],
                        summary(EE_3id)$coefficients[4,1], summary(EE_3id)$coefficients[4,2],
                        summary(EE_3id)$coefficients[5,1], summary(EE_3id)$coefficients[5,2],
                        summary(EE_3id)$coefficients[6,1], summary(EE_3id)$coefficients[6,2],
                        logLik(EE_3id), with(summary(EE_3id), 1-deviance/null.deviance)),digits=2)

TableA19


##++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# C.4 Separated Logistic Regressions: Natives and Non-natives
##++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#===================================================================================
# Estimations (for Table A21)
#===================================================================================

#---------------------------
# Germany (DE)
#---------------------------
# subsample foreigner
DE_3f <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
             data=DE[DE$foreigner=="1", ], family="binomial")
summary(DE_3f)
# subsample non-foreigner
DE_3nf <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben, 
              data=DE[DE$foreigner=="0", ], family="binomial")
summary(DE_3nf)

#---------------------------
# France (FR)
#---------------------------
# subsample foreigner
FR_3f <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
             data=FR[FR$foreigner=="1", ], family="binomial")
summary(FR_3f)
# subsample non-foreigner
FR_3nf <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
              data=FR[FR$foreigner=="0", ], family="binomial")
summary(FR_3nf)

#---------------------------
# Italy (IT)
#---------------------------
# subsample foreigner
IT_3f <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
             data=IT[IT$foreigner=="1", ], family="binomial")
summary(IT_3f)
# subsample non-foreigner
IT_3nf <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
              data=IT[IT$foreigner=="0", ], family="binomial")
summary(IT_3nf)

#---------------------------
# Spain (ES)
#---------------------------
# subsample foreigner
ES_3f <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
             data=ES[ES$foreigner=="1", ], family="binomial")
summary(ES_3f)
# subsample non-foreigner
ES_3nf <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
              data=ES[ES$foreigner=="0", ], family="binomial")
summary(ES_3nf)

#---------------------------
# Poland (PL)
#---------------------------
# subsample foreigner
PL_3f <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
             data=PL[PL$foreigner=="1", ], family="binomial")
summary(PL_3f)
# subsample non-foreigner
PL_3nf <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
              data=PL[PL$foreigner=="0", ], family="binomial")
summary(PL_3nf)

#---------------------------
# Romania (RO)
#---------------------------
# subsample foreigner
RO_3f <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
             data=RO[RO$foreigner=="1", ], family="binomial")
summary(RO_3f)
# subsample non-foreigner
RO_3nf <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
              data=RO[RO$foreigner=="0", ], family="binomial")
summary(RO_3nf)

#---------------------------
# Netherlands (NL)
#---------------------------
# subsample foreigner
NL_3f <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
             data=NL[NL$foreigner=="1", ], family="binomial")
summary(NL_3f)
# subsample non-foreigner
NL_3nf <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
              data=NL[NL$foreigner=="0", ], family="binomial")
summary(NL_3nf)

#---------------------------
# Sweden (SE)
#---------------------------
# subsample foreigner
SE_3f <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
             data=SE[SE$foreigner=="1", ], family="binomial")
summary(SE_3f)
# subsample non-foreigner
SE_3nf <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
              data=SE[SE$foreigner=="0", ], family="binomial")
summary(SE_3nf)

#---------------------------
# Hungary (HU)
#---------------------------
# subsample foreigner
HU_3f <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
             data=HU[HU$foreigner=="1", ], family="binomial")
summary(HU_3f)
# subsample non-foreigner
HU_3nf <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
              data=HU[HU$foreigner=="0", ], family="binomial")
summary(HU_3nf)

#---------------------------
# Austria (AT)
#---------------------------
# subsample foreigner
AT_3f <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
             data=AT[AT$foreigner=="1", ], family="binomial")
summary(AT_3f)
# subsample non-foreigner
AT_3nf <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
              data=AT[AT$foreigner=="0", ], family="binomial")
summary(AT_3nf)

#---------------------------
# Bulgaria (BG)
#---------------------------
# subsample foreigner
BG_3f <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
             data=BG[BG$foreigner=="1", ], family="binomial")
summary(BG_3f)
# subsample non-foreigner
BG_3nf <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
              data=BG[BG$foreigner=="0", ], family="binomial")
summary(BG_3nf)

#---------------------------
# Slovakia (SK)
#---------------------------
# subsample foreigner
SK_3f <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
             data=SK[SK$foreigner=="1", ], family="binomial")
summary(SK_3f)
# subsample non-foreigner
SK_3nf <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
              data=SK[SK$foreigner=="0", ], family="binomial")
summary(SK_3nf)

#---------------------------
# Latvia (LV)
#---------------------------
# subsample foreigner
LV_3f <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
             data=LV[LV$foreigner=="1", ], family="binomial")
summary(LV_3f)
# subsample non-foreigner
LV_3nf <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
              data=LV[LV$foreigner=="0", ], family="binomial")
summary(LV_3nf)

#---------------------------
# Estonia (EE)
#---------------------------
# subsample foreigner (not truth worthy due to no variation in the variable across the response; not reported)
EE_3f <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
             data=EE[EE$foreigner=="1", ], family="binomial")
summary(EE_3f)

# subsample non-foreigner
EE_3nf <- glm(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben,
              data=EE[EE$foreigner=="0", ], family="binomial")
summary(EE_3nf)


#------------------------------------------------------------------------------------
# Table A20: Estimates on Ideological and Identity Variables: Natives and Non-natives
#------------------------------------------------------------------------------------

TableA20 <- matrix(NA, 42, ncol = 7)
colnames(TableA20) <- c("beta1", "beta2", "beta3", "beta4", "Log-Lik.", "pse. R^2","Obs.")
rownames(TableA20) <- c("Germany: Complete Model", "Germany: Non-Native Model", "Germany: Native Model", 
                        "France: Complete Model", "France: Non-Native Model", "France: Native Model",
                        "Italy: Complete Model", "Italy: Non-Native Model", "Italy: Native Model",
                        "Spain: Complete Model", "Spain: Non-Native Model", "Spain: Native Model",
                        "Poland: Complete Model", "Poland: Non-Native Model", "Poland: Native Model",
                        "Romania: Complete Model", "Romania: Non-Native Model", "Romania: Native Model",
                        "Netherlands: Complete Model", "Netherlands: Non-Native Model", "Netherlands: Native Model",
                        "Sweden: Complete Model", "Sweden: Non-Native Model", "Sweden: Native Model",
                        "Hungary: Complete Model", "Hungary: Non-Native Model", "Hungary: Native Model",
                        "Austria: Complete Model", "Austria: Non-Native Model", "Austria: Native Model",
                        "Bulgaria: Complete Model", "Bulgaria: Non-Native Model", "Bulgaria: Native Model",
                        "Slovakia: Complete Model", "Slovakia: Non-Native Model", "Slovakia: Native Model",
                        "Latvia: Complete Model", "Latvia: Non-Native Model", "Latvia: Native Model",
                        "Estonia: Complete Model", "Estonia: Non-Native Model", "Estonia: Native Model")
# DE
TableA20[1,]  <- c(round(summary(DE_3)$coefficients[2,1], digits=2),
                   round(summary(DE_3)$coefficients[3,1], digits=2), 
                   round(summary(DE_3)$coefficients[4,1], digits=2), 
                   round(summary(DE_3)$coefficients[5,1], digits=2), 
                   round(logLik(DE_3), digits=2), round(with(summary(DE_3), 1 - deviance/null.deviance), digits=2),nobs(DE_3))

TableA20[2,]  <- c(round(summary(DE_3f)$coefficients[2,1], digits=2), 
                   round(summary(DE_3f)$coefficients[3,1], digits=2), 
                   round(summary(DE_3f)$coefficients[4,1], digits=2),
                   round(summary(DE_3f)$coefficients[5,1], digits=2), 
                   round(logLik(DE_3f), digits=2), round(with(summary(DE_3f), 1 - deviance/null.deviance), digits=2),nobs(DE_3f))

TableA20[3,]  <- c(round(summary(DE_3nf)$coefficients[2,1], digits=2), 
                   round(summary(DE_3nf)$coefficients[3,1], digits=2), 
                   round(summary(DE_3nf)$coefficients[4,1], digits=2), 
                   round(summary(DE_3nf)$coefficients[5,1], digits=2), 
                   round(logLik(DE_3nf), digits=2), round(with(summary(DE_3nf), 1 - deviance/null.deviance), digits=2),nobs(DE_3nf))
# FR
TableA20[4,]  <- c(round(summary(FR_3)$coefficients[2,1], digits=2), 
                   round(summary(FR_3)$coefficients[3,1], digits=2),
                   round(summary(FR_3)$coefficients[4,1], digits=2), 
                   round(summary(FR_3)$coefficients[5,1], digits=2), 
                   round(logLik(FR_3), digits=2), round(with(summary(FR_3), 1 - deviance/null.deviance), digits=2),
                   nobs(FR_3))

TableA20[5,]  <- c(round(summary(FR_3f)$coefficients[2,1], digits=2), 
                   round(summary(FR_3f)$coefficients[3,1], digits=2), 
                   round(summary(FR_3f)$coefficients[4,1], digits=2), 
                   round(summary(FR_3f)$coefficients[5,1], digits=2), 
                   round(logLik(FR_3f), digits=2), round(with(summary(FR_3f), 1 - deviance/null.deviance), digits=2),
                   nobs(FR_3f))

TableA20[6,]  <- c(round(summary(FR_3nf)$coefficients[2,1], digits=2), 
                   round(summary(FR_3nf)$coefficients[3,1], digits=2), 
                   round(summary(FR_3nf)$coefficients[4,1], digits=2), 
                   round(summary(FR_3nf)$coefficients[5,1], digits=2), 
                   round(logLik(FR_3nf), digits=2), round(with(summary(FR_3nf), 1 - deviance/null.deviance), digits=2),
                   nobs(FR_3nf))
# IT
TableA20[7,]  <- c(round(summary(IT_3)$coefficients[2,1], digits=2),
                   round(summary(IT_3)$coefficients[3,1], digits=2),
                   round(summary(IT_3)$coefficients[4,1], digits=2), 
                   round(summary(IT_3)$coefficients[5,1], digits=2), 
                   round(logLik(IT_3), digits=2), round(with(summary(IT_3), 1 - deviance/null.deviance), digits=2),
                   nobs(IT_3))

TableA20[8,]  <- c(round(summary(IT_3f)$coefficients[2,1], digits=2), 
                   round(summary(IT_3f)$coefficients[3,1], digits=2), 
                   round(summary(IT_3f)$coefficients[4,1], digits=2), 
                   round(summary(IT_3f)$coefficients[5,1], digits=2), 
                   round(logLik(IT_3f), digits=2), round(with(summary(IT_3f), 1 - deviance/null.deviance), digits=2),
                   nobs(IT_3f))

TableA20[9,]  <- c(round(summary(IT_3nf)$coefficients[2,1], digits=2), 
                   round(summary(IT_3nf)$coefficients[3,1], digits=2), 
                   round(summary(IT_3nf)$coefficients[4,1], digits=2),
                   round(summary(IT_3nf)$coefficients[5,1], digits=2), 
                   round(logLik(IT_3nf), digits=2), round(with(summary(IT_3nf), 1 - deviance/null.deviance), digits=2),
                   nobs(IT_3nf))
# ES
TableA20[10,] <- c(round(summary(ES_3)$coefficients[2,1], digits=2), 
                   round(summary(ES_3)$coefficients[3,1], digits=2), 
                   round(summary(ES_3)$coefficients[4,1], digits=2), 
                   round(summary(ES_3)$coefficients[5,1], digits=2), 
                   round(logLik(ES_3), digits=2), round(with(summary(ES_3), 1 - deviance/null.deviance), digits=2),
                   nobs(ES_3))

TableA20[11,] <- c(round(summary(ES_3f)$coefficients[2,1], digits=2), 
                   round(summary(ES_3f)$coefficients[3,1], digits=2), 
                   round(summary(ES_3f)$coefficients[4,1], digits=2), 
                   round(summary(ES_3f)$coefficients[5,1], digits=2), 
                   round(logLik(ES_3f), digits=2), round(with(summary(ES_3f), 1 - deviance/null.deviance), digits=2),
                   nobs(ES_3f))

TableA20[12,] <- c(round(summary(ES_3nf)$coefficients[2,1], digits=2), 
                   round(summary(ES_3nf)$coefficients[3,1], digits=2), 
                   round(summary(ES_3nf)$coefficients[4,1], digits=2), 
                   round(summary(ES_3nf)$coefficients[5,1], digits=2), 
                   round(logLik(ES_3nf), digits=2), round(with(summary(ES_3nf), 1 - deviance/null.deviance), digits=2),
                   nobs(ES_3nf))
# PL
TableA20[13,]  <-c(round(summary(PL_3)$coefficients[2,1], digits=2), 
                   round(summary(PL_3)$coefficients[3,1], digits=2), 
                   round(summary(PL_3)$coefficients[4,1], digits=2), 
                   round(summary(PL_3)$coefficients[5,1], digits=2),
                   round(logLik(PL_3), digits=2), round(with(summary(PL_3), 1  -deviance/null.deviance), digits=2),nobs(PL_3))

TableA20[14,]  <-c(round(summary(PL_3f)$coefficients[2,1], digits=2), 
                   round(summary(PL_3f)$coefficients[3,1], digits=2), 
                   round(summary(PL_3f)$coefficients[4,1], digits=2),
                   round(summary(PL_3f)$coefficients[5,1], digits=2), 
                   round(logLik(PL_3f), digits=2), round(with(summary(PL_3f), 1  -deviance/null.deviance), digits=2),nobs(PL_3f))

TableA20[15,]  <-c(round(summary(PL_3nf)$coefficients[2,1], digits=2), 
                   round(summary(PL_3nf)$coefficients[3,1], digits=2), 
                   round(summary(PL_3nf)$coefficients[4,1], digits=2), 
                   round(summary(PL_3nf)$coefficients[5,1], digits=2), 
                   round(logLik(PL_3nf), digits=2), round(with(summary(PL_3nf), 1  -deviance/null.deviance), digits=2),nobs(PL_3nf))
# RO
TableA20[16,]  <-c(round(summary(RO_3)$coefficients[2,1], digits=2), 
                   round(summary(RO_3)$coefficients[3,1], digits=2), 
                   round(summary(RO_3)$coefficients[4,1], digits=2), 
                   round(summary(RO_3)$coefficients[5,1], digits=2), 
                   round(logLik(RO_3), digits=2), round(with(summary(RO_3), 1  -deviance/null.deviance), digits=2),nobs(RO_3))

TableA20[17,]  <-c(round(summary(RO_3f)$coefficients[2,1], digits=2), 
                   round(summary(RO_3f)$coefficients[3,1], digits=2), 
                   round(summary(RO_3f)$coefficients[4,1], digits=2), 
                   round(summary(RO_3f)$coefficients[5,1], digits=2), 
                   round(logLik(RO_3f), digits=2), round(with(summary(RO_3f), 1  -deviance/null.deviance), digits=2),nobs(RO_3f))

TableA20[18,]  <-c(round(summary(RO_3nf)$coefficients[2,1], digits=2),
                   round(summary(RO_3nf)$coefficients[3,1], digits=2), 
                   round(summary(RO_3nf)$coefficients[4,1], digits=2),
                   round(summary(RO_3nf)$coefficients[5,1], digits=2), 
                   round(logLik(RO_3nf), digits=2), round(with(summary(RO_3nf), 1  -deviance/null.deviance), digits=2),nobs(RO_3nf))
# NL
TableA20[19,]  <-c(round(summary(NL_3)$coefficients[2,1], digits=2), 
                   round(summary(NL_3)$coefficients[3,1], digits=2), 
                   round(summary(NL_3)$coefficients[4,1], digits=2), 
                   round(summary(NL_3)$coefficients[5,1], digits=2), 
                   round(logLik(NL_3), digits=2), round(with(summary(NL_3), 1  -deviance/null.deviance), digits=2),nobs(NL_3))

TableA20[20,]  <-c(round(summary(NL_3f)$coefficients[2,1], digits=2), 
                   round(summary(NL_3f)$coefficients[3,1], digits=2), 
                   round(summary(NL_3f)$coefficients[4,1], digits=2), 
                   round(summary(NL_3f)$coefficients[5,1], digits=2), 
                   round(logLik(NL_3f), digits=2), round(with(summary(NL_3f), 1  -deviance/null.deviance), digits=2),nobs(NL_3f))

TableA20[21,]  <-c(round(summary(NL_3nf)$coefficients[2,1], digits=2), 
                   round(summary(NL_3nf)$coefficients[3,1], digits=2),
                   round(summary(NL_3nf)$coefficients[4,1], digits=2), 
                   round(summary(NL_3nf)$coefficients[5,1], digits=2), 
                   round(logLik(NL_3nf), digits=2), round(with(summary(NL_3nf), 1  -deviance/null.deviance), digits=2),nobs(NL_3nf))
# SE
TableA20[22,]  <-c(round(summary(SE_3)$coefficients[2,1], digits=2), 
                   round(summary(SE_3)$coefficients[3,1], digits=2), 
                   round(summary(SE_3)$coefficients[4,1], digits=2), 
                   round(summary(SE_3)$coefficients[5,1], digits=2), 
                   round(logLik(SE_3), digits=2), round(with(summary(SE_3), 1  -deviance/null.deviance), digits=2),nobs(SE_3))

TableA20[23,]  <-c(round(summary(SE_3f)$coefficients[2,1], digits=2), 
                   round(summary(SE_3f)$coefficients[3,1], digits=2), 
                   round(summary(SE_3f)$coefficients[4,1], digits=2), 
                   round(summary(SE_3f)$coefficients[5,1], digits=2), 
                   round(logLik(SE_3f), digits=2), round(with(summary(SE_3f), 1  -deviance/null.deviance), digits=2),nobs(SE_3f))

TableA20[24,]  <-c(round(summary(SE_3nf)$coefficients[2,1], digits=2), 
                   round(summary(SE_3nf)$coefficients[3,1], digits=2), 
                   round(summary(SE_3nf)$coefficients[4,1], digits=2), 
                   round(summary(SE_3nf)$coefficients[5,1], digits=2), 
                   round(logLik(SE_3nf), digits=2), round(with(summary(SE_3nf), 1  -deviance/null.deviance), digits=2),nobs(SE_3nf))
# HU
TableA20[25,]  <-c(round(summary(HU_3)$coefficients[2,1], digits=2), 
                   round(summary(HU_3)$coefficients[3,1], digits=2), 
                   round(summary(HU_3)$coefficients[4,1], digits=2), 
                   round(summary(HU_3)$coefficients[5,1], digits=2), 
                   round(logLik(HU_3), digits=2), round(with(summary(HU_3), 1  -deviance/null.deviance), digits=2),nobs(HU_3))

TableA20[26,]  <-c(round(summary(HU_3f)$coefficients[2,1], digits=2), 
                   round(summary(HU_3f)$coefficients[3,1], digits=2), 
                   round(summary(HU_3f)$coefficients[4,1], digits=2), 
                   round(summary(HU_3f)$coefficients[5,1], digits=2),
                   round(logLik(HU_3f), digits=2), round(with(summary(HU_3f), 1  -deviance/null.deviance), digits=2),nobs(HU_3f))

TableA20[27,]  <-c(round(summary(HU_3nf)$coefficients[2,1], digits=2), 
                   round(summary(HU_3nf)$coefficients[3,1], digits=2), 
                   round(summary(HU_3nf)$coefficients[4,1], digits=2), 
                   round(summary(HU_3nf)$coefficients[5,1], digits=2), 
                   round(logLik(HU_3nf), digits=2), round(with(summary(HU_3nf), 1  -deviance/null.deviance), digits=2),nobs(HU_3nf))
# AT
TableA20[28,]  <-c(round(summary(AT_3)$coefficients[2,1], digits=2), 
                   round(summary(AT_3)$coefficients[3,1], digits=2), 
                   round(summary(AT_3)$coefficients[4,1], digits=2), 
                   round(summary(AT_3)$coefficients[5,1], digits=2), 
                   round(logLik(AT_3), digits=2), round(with(summary(AT_3), 1  -deviance/null.deviance), digits=2),nobs(AT_3))

TableA20[29,]  <-c(round(summary(AT_3f)$coefficients[2,1], digits=2), 
                   round(summary(AT_3f)$coefficients[3,1], digits=2), 
                   round(summary(AT_3f)$coefficients[4,1], digits=2), 
                   round(summary(AT_3f)$coefficients[5,1], digits=2), 
                   round(logLik(AT_3f), digits=2), round(with(summary(AT_3f), 1  -deviance/null.deviance), digits=2),nobs(AT_3f))

TableA20[30,]  <-c(round(summary(AT_3nf)$coefficients[2,1], digits=2), 
                   round(summary(AT_3nf)$coefficients[3,1], digits=2), 
                   round(summary(AT_3nf)$coefficients[4,1], digits=2), 
                   round(summary(AT_3nf)$coefficients[5,1], digits=2), 
                   round(logLik(AT_3nf), digits=2), round(with(summary(AT_3nf), 1  -deviance/null.deviance), digits=2),nobs(AT_3nf))
# BG
TableA20[31,]  <-c(round(summary(BG_3)$coefficients[2,1], digits=2), 
                   round(summary(BG_3)$coefficients[3,1], digits=2),
                   round(summary(BG_3)$coefficients[4,1], digits=2), 
                   round(summary(BG_3)$coefficients[5,1], digits=2), 
                   round(logLik(BG_3), digits=2), round(with(summary(BG_3), 1  -deviance/null.deviance), digits=2),nobs(BG_3))

TableA20[32,]  <-c(round(summary(BG_3f)$coefficients[2,1], digits=2), 
                   round(summary(BG_3f)$coefficients[3,1], digits=2), 
                   round(summary(BG_3f)$coefficients[4,1], digits=2), 
                   round(summary(BG_3f)$coefficients[5,1], digits=2), 
                   round(logLik(BG_3f), digits=2), round(with(summary(BG_3f), 1  -deviance/null.deviance), digits=2),nobs(BG_3f))

TableA20[33,]  <-c(round(summary(BG_3nf)$coefficients[2,1], digits=2), 
                   round(summary(BG_3nf)$coefficients[3,1], digits=2), 
                   round(summary(BG_3nf)$coefficients[4,1], digits=2), 
                   round(summary(BG_3nf)$coefficients[5,1], digits=2), 
                   round(logLik(BG_3nf), digits=2), round(with(summary(BG_3nf), 1  -deviance/null.deviance), digits=2),nobs(BG_3nf))
# SK
TableA20[34,]  <-c(round(summary(SK_3)$coefficients[2,1], digits=2), 
                   round(summary(SK_3)$coefficients[3,1], digits=2), 
                   round(summary(SK_3)$coefficients[4,1], digits=2), 
                   round(summary(SK_3)$coefficients[5,1], digits=2), 
                   round(logLik(SK_3), digits=2), round(with(summary(SK_3), 1  -deviance/null.deviance), digits=2),nobs(SK_3))

TableA20[35,]  <-c(round(summary(SK_3f)$coefficients[2,1], digits=2), 
                   round(summary(SK_3f)$coefficients[3,1], digits=2), 
                   round(summary(SK_3f)$coefficients[4,1], digits=2), 
                   round(summary(SK_3f)$coefficients[5,1], digits=2), 
                   round(logLik(SK_3f), digits=2), round(with(summary(SK_3f), 1  -deviance/null.deviance), digits=2),nobs(SK_3f))

TableA20[36,]  <-c(round(summary(SK_3nf)$coefficients[2,1], digits=2), 
                   round(summary(SK_3nf)$coefficients[3,1], digits=2), 
                   round(summary(SK_3nf)$coefficients[4,1], digits=2), 
                   round(summary(SK_3nf)$coefficients[5,1], digits=2),
                   round(logLik(SK_3nf), digits=2), round(with(summary(SK_3nf), 1  -deviance/null.deviance), digits=2),nobs(SK_3nf))
# LV
TableA20[37,]  <-c(round(summary(LV_3)$coefficients[2,1], digits=2), 
                   round(summary(LV_3)$coefficients[3,1], digits=2), 
                   round(summary(LV_3)$coefficients[4,1], digits=2), 
                   round(summary(LV_3)$coefficients[5,1], digits=2), 
                   round(logLik(LV_3), digits=2), round(with(summary(LV_3), 1  -deviance/null.deviance), digits=2),nobs(LV_3))

TableA20[38,]  <-c(round(summary(LV_3f)$coefficients[2,1], digits=2), 
                   round(summary(LV_3f)$coefficients[3,1], digits=2), 
                   round(summary(LV_3f)$coefficients[4,1], digits=2), 
                   round(summary(LV_3f)$coefficients[5,1], digits=2), 
                   round(logLik(LV_3f), digits=2), round(with(summary(LV_3f), 1  -deviance/null.deviance), digits=2),nobs(LV_3f))

TableA20[39,]  <-c(round(summary(LV_3nf)$coefficients[2,1], digits=2), 
                   round(summary(LV_3nf)$coefficients[3,1], digits=2), 
                   round(summary(LV_3nf)$coefficients[4,1], digits=2), 
                   round(summary(LV_3nf)$coefficients[5,1], digits=2), 
                   round(logLik(LV_3nf), digits=2), round(with(summary(LV_3nf), 1  -deviance/null.deviance), digits=2),nobs(LV_3nf))
# EE
TableA20[40,]  <-c(round(summary(EE_3)$coefficients[2,1], digits=2), 
                   round(summary(EE_3)$coefficients[3,1], digits=2), 
                   round(summary(EE_3)$coefficients[4,1], digits=2), 
                   round(summary(EE_3)$coefficients[5,1], digits=2), 
                   round(logLik(EE_3), digits=2), round(with(summary(EE_3), 1 -deviance/null.deviance), digits=2),nobs(EE_3)
)

TableA20[41,]  <-c(
  round(summary(EE_3f)$coefficients[2,1], digits=2), 
  round(summary(EE_3f)$coefficients[3,1], digits=2), 
  round(summary(EE_3f)$coefficients[4,1], digits=2), 
  round(summary(EE_3f)$coefficients[5,1], digits=2), 
  round(logLik(EE_3f), digits=2), round(with(summary(EE_3f), 1 -deviance/null.deviance), digits=2),nobs(EE_3f)
)

TableA20[42,]  <-c(
  round(summary(EE_3nf)$coefficients[2,1], digits=2),
  round(summary(EE_3nf)$coefficients[3,1], digits=2), 
  round(summary(EE_3nf)$coefficients[4,1], digits=2),
  round(summary(EE_3nf)$coefficients[5,1], digits=2), 
  round(logLik(EE_3nf), digits=2), round(with(summary(EE_3nf), 1 -deviance/null.deviance), digits=2),nobs(EE_3nf)
)


TableA20


##++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# C.5 Pooled and Multilevel Models 
##++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#===================================================================================
# Table A21: Binary Logistic Regression Model: Pooled 
#===================================================================================

TableA21 <- matrix(NA, 15, ncol = 3)
colnames(TableA21) <- c("coef.", "s.e.", "p-value" )
rownames(TableA21) <- c( "Intercept", "Immigration", "Redistribution", "National-European", "National-Regional",
                         "Female", "Age", "Population Size: small city", "Population Size: large city",
                         "Education: high school", "Education: university","Economic Situation: satisfied",
                         "Income: medium", "Income: high", "Benefit EU Funds: yes")
TableA21[,1]   <- round(summary(FS_P)$coefficients[,1], digits=2)
TableA21[,2]   <- round(summary(FS_P)$coefficients[,2], digits=2)
TableA21[,3]   <- round(summary(FS_P)$coefficients[,4], digits=2)
TableA21

#===================================================================================
# Estimation Random Intercept (for Table A22)
#===================================================================================

FS_mlri <- glmer(eu_support ~ pro_imms + pro_reds + nat_eus + nat_regs + female + aged + pop + educ + eco + income + eu_ben + (1 | country), data=FS, family="binomial")
summary(FS_mlri)

#===================================================================================
# Table A22: Binary Logistic Regression Model: Random Intercept  
#===================================================================================

r_intercept <- as.data.frame(VarCorr(FS_mlri)) #variance and standard deviation random intercept

TableA22a <- matrix(NA, 15, ncol = 3)
colnames(TableA22a) <- c("coef.", "s.e.", "p-value" )
rownames(TableA22a) <- c( "Intercept", "Immigration", "Redistribution", "National-European", "National-Regional",
                          "Female", "Age", "Population Size: small city", "Population Size: large city",
                          "Education: high school", "Education: university","Economic Situation: satisfied",
                          "Income: medium", "Income: high", "Benefit EU Funds: yes")


TableA22a[,1]   <- round(summary(FS_mlri)$coefficients[,1], digits=2)
TableA22a[,2]   <- round(summary(FS_mlri)$coefficients[,2], digits=2)
TableA22a[,3]   <- round(summary(FS_mlri)$coefficients[,4], digits=2)
TableA22a

TableA22b <- matrix(NA, 1, ncol = 2)
colnames(TableA22b) <- c("var.", "st.dev.")
rownames(TableA22b) <- c( "Intercept")

TableA22b[1,]   <- round(c(r_intercept[,4], r_intercept[,5]),digits=2)

TableA22b


#===================================================================================
# Model Comparisons (for Table A23)
#===================================================================================

# Loglikelihood Pooled vs Random Intercepts

ll0 <- as.numeric(logLik(FS_P))   # logLik null model (pooled)
ll1 <- as.numeric(logLik(FS_mlri))    # logLik full model (random intercepts)
LR <- -2 * (ll0 - ll1)
p_mixed <- 0.5 * pchisq(LR, df = 1, lower.tail = FALSE)
LR; p_mixed

# AIC
AIC(FS_P)
AIC(FS_mlri)

# BIC
BIC(FS_P)
BIC(FS_mlri)


#===================================================================================
# Table A23: Explaining Support for EU Membership: Model Comparisons 
#===================================================================================

TableA23 <- matrix(NA, 2, ncol = 3)
colnames(TableA23) <- c("Log-Lik", "AIC", "BIC")
rownames(TableA23) <- c("Pooled", "Random Intercepts")

TableA23[1,]<-round(c(logLik(FS_P), AIC(FS_P), BIC(FS_P)),digits=2)
TableA23[2,]<-round(c(logLik(FS_mlri), AIC(FS_mlri), BIC(FS_mlri)),digits=2)

TableA23

#-------------------------------------
# D Robustness Checks
#-------------------------------------

##++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# D.1 Comparing the Binary and the Cumulative Logistic Regression Models
##++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


#------------------------------------------------------------------------------------
# Table A24: Robustness Check: Summary Table. Binary and Cumulative Logistic Models
#------------------------------------------------------------------------------------

TableA24 <- matrix(NA, 28, ncol = 7)
colnames(TableA24) <- c("Obs.", "beta1", "beta2", "beta3", "beta4", "Log-Lik", "pse. R^2")
rownames(TableA24) <- c("Germany: binary", "Germany: cumulative", "France: binary", "France: cumulative", "Italy: binary", "Italy: cumulative", 
                        "Spain: binary", "Spain: cumulative", "Poland: binary", "Poland: cumulative", "Romania: binary", "Romania: cumulative", 
                        "Netherlands: binary", "Netherlands: cumulative", "Sweden: binary",  "Sweden: cumulative", "Hungary: binary", "Hungary: cumulative",
                        "Austria: binary", "Austria: cumulative", "Bulgaria: binary", "Bulgaria: cumulative", "Slovakia: binary", "Slovakia: cumulative",
                        "Latvia: binary", "Latvia: cumulative", "Estonia: binary", "Estonia: cumulative")

# DE
TableA24[1,] <-round(c(nrow(DE),summary(DE_3)$coefficients[2,1], 
                       summary(DE_3)$coefficients[3,1], 
                       summary(DE_3)$coefficients[4,1], 
                       summary(DE_3)$coefficients[5,1], 
                       logLik(DE_3), with(summary(DE_3), 1-deviance/null.deviance)),digits=2)

TableA24[2,] <-round(c(nrow(DE3),coef(DE_3_cum)["pro_imms"],
                       coef(DE_3_cum)["pro_reds"],
                       coef(DE_3_cum)["nat_eus"], 
                       coef(DE_3_cum)["nat_regs"],
                       logLik(DE_3_cum), 1 - (logLik(DE_3_cum) / as.numeric(logLik(update(DE_3_cum, . ~ 1))))), digits = 2)
# FR
TableA24[3,] <-round(c(nrow(FR),summary(FR_3)$coefficients[2,1], 
                       summary(FR_3)$coefficients[3,1], 
                       summary(FR_3)$coefficients[4,1], 
                       summary(FR_3)$coefficients[5,1], 
                       logLik(FR_3), with(summary(FR_3), 1-deviance/null.deviance)),digits=2)

TableA24[4,] <-round(c(nrow(FR3),coef(FR_3_cum)["pro_imms"],
                       coef(FR_3_cum)["pro_reds"],
                       coef(FR_3_cum)["nat_eus"], 
                       coef(FR_3_cum)["nat_regs"],
                       logLik(FR_3_cum), 1 - (logLik(FR_3_cum) / as.numeric(logLik(update(FR_3_cum, . ~ 1))))), digits = 2)
# IT
TableA24[5,] <-round(c(nrow(IT),summary(IT_3)$coefficients[2,1], 
                       summary(IT_3)$coefficients[3,1], 
                       summary(IT_3)$coefficients[4,1], 
                       summary(IT_3)$coefficients[5,1], 
                       logLik(IT_3), with(summary(IT_3), 1-deviance/null.deviance)),digits=2)

TableA24[6,] <-round(c(nrow(IT3),coef(IT_3_cum)["pro_imms"],
                       coef(IT_3_cum)["pro_reds"],
                       coef(IT_3_cum)["nat_eus"], 
                       coef(IT_3_cum)["nat_regs"],
                       logLik(IT_3_cum), 1 - (logLik(IT_3_cum) / as.numeric(logLik(update(IT_3_cum, . ~ 1))))), digits = 2)
# ES
TableA24[7,] <-round(c(nrow(ES),summary(ES_3)$coefficients[2,1], 
                       summary(ES_3)$coefficients[3,1], 
                       summary(ES_3)$coefficients[4,1], 
                       summary(ES_3)$coefficients[5,1], 
                       logLik(ES_3), with(summary(ES_3), 1-deviance/null.deviance)),digits=2)

TableA24[8,] <-round(c(nrow(ES3),coef(ES_3_cum)["pro_imms"],
                       coef(ES_3_cum)["pro_reds"],
                       coef(ES_3_cum)["nat_eus"], 
                       coef(ES_3_cum)["nat_regs"],
                       logLik(ES_3_cum), 1 - (logLik(ES_3_cum) / as.numeric(logLik(update(ES_3_cum, . ~ 1))))), digits = 2)
# PL
TableA24[9,] <-round(c(nrow(PL),summary(PL_3)$coefficients[2,1], 
                       summary(PL_3)$coefficients[3,1], 
                       summary(PL_3)$coefficients[4,1], 
                       summary(PL_3)$coefficients[5,1], 
                       logLik(PL_3), with(summary(PL_3), 1-deviance/null.deviance)),digits=2)

TableA24[10,] <-round(c(nrow(PL3),coef(PL_3_cum)["pro_imms"],
                        coef(PL_3_cum)["pro_reds"],
                        coef(PL_3_cum)["nat_eus"], 
                        coef(PL_3_cum)["nat_regs"],
                        logLik(PL_3_cum), 1 - (logLik(PL_3_cum) / as.numeric(logLik(update(PL_3_cum, . ~ 1))))), digits = 2)
# RO
TableA24[11,] <-round(c(nrow(RO),summary(RO_3)$coefficients[2,1], 
                        summary(RO_3)$coefficients[3,1], 
                        summary(RO_3)$coefficients[4,1], 
                        summary(RO_3)$coefficients[5,1], 
                        logLik(RO_3), with(summary(RO_3), 1-deviance/null.deviance)),digits=2)

TableA24[12,] <-round(c(nrow(RO3),coef(RO_3_cum)["pro_imms"],
                        coef(RO_3_cum)["pro_reds"],
                        coef(RO_3_cum)["nat_eus"], 
                        coef(RO_3_cum)["nat_regs"],
                        logLik(RO_3_cum), 1 - (logLik(RO_3_cum) / as.numeric(logLik(update(RO_3_cum, . ~ 1))))), digits = 2)
# NL
TableA24[13,] <-round(c(nrow(NL),summary(NL_3)$coefficients[2,1], 
                        summary(NL_3)$coefficients[3,1], 
                        summary(NL_3)$coefficients[4,1], 
                        summary(NL_3)$coefficients[5,1], 
                        logLik(NL_3), with(summary(NL_3), 1-deviance/null.deviance)),digits=2)

TableA24[14,] <-round(c(nrow(NL3),coef(NL_3_cum)["pro_imms"],
                        coef(NL_3_cum)["pro_reds"],
                        coef(NL_3_cum)["nat_eus"], 
                        coef(NL_3_cum)["nat_regs"],
                        logLik(NL_3_cum), 1 - (logLik(NL_3_cum) / as.numeric(logLik(update(NL_3_cum, . ~ 1))))), digits = 2)
# SE
TableA24[15,] <-round(c(nrow(SE),summary(SE_3)$coefficients[2,1],
                        summary(SE_3)$coefficients[3,1], 
                        summary(SE_3)$coefficients[4,1], 
                        summary(SE_3)$coefficients[5,1], 
                        logLik(SE_3), with(summary(SE_3), 1-deviance/null.deviance)),digits=2)

TableA24[16,] <-round(c(nrow(SE3),coef(SE_3_cum)["pro_imms"],
                        coef(SE_3_cum)["pro_reds"],
                        coef(SE_3_cum)["nat_eus"], 
                        coef(SE_3_cum)["nat_regs"],
                        logLik(SE_3_cum), 1 - (logLik(SE_3_cum) / as.numeric(logLik(update(SE_3_cum, . ~ 1))))), digits = 2)
# HU
TableA24[17,] <-round(c(nrow(HU),summary(HU_3)$coefficients[2,1], 
                        summary(HU_3)$coefficients[3,1], 
                        summary(HU_3)$coefficients[4,1], 
                        summary(HU_3)$coefficients[5,1], 
                        logLik(HU_3), with(summary(HU_3), 1-deviance/null.deviance)),digits=2)

TableA24[18,] <-round(c(nrow(HU3),coef(HU_3_cum)["pro_imms"],
                        coef(HU_3_cum)["pro_reds"],
                        coef(HU_3_cum)["nat_eus"], 
                        coef(HU_3_cum)["nat_regs"],
                        logLik(HU_3_cum), 1 - (logLik(HU_3_cum) / as.numeric(logLik(update(HU_3_cum, . ~ 1))))), digits = 2)
# AT
TableA24[19,] <-round(c(nrow(AT),summary(AT_3)$coefficients[2,1], 
                        summary(AT_3)$coefficients[3,1], 
                        summary(AT_3)$coefficients[4,1], 
                        summary(AT_3)$coefficients[5,1], 
                        logLik(AT_3), with(summary(AT_3), 1-deviance/null.deviance)),digits=2)

TableA24[20,] <-round(c(nrow(AT3),coef(AT_3_cum)["pro_imms"],
                        coef(AT_3_cum)["pro_reds"],
                        coef(AT_3_cum)["nat_eus"], 
                        coef(AT_3_cum)["nat_regs"],
                        logLik(AT_3_cum), 1 - (logLik(AT_3_cum) / as.numeric(logLik(update(AT_3_cum, . ~ 1))))), digits = 2)
# BG
TableA24[21,] <-round(c(nrow(BG),summary(BG_3)$coefficients[2,1], 
                        summary(BG_3)$coefficients[3,1], 
                        summary(BG_3)$coefficients[4,1], 
                        summary(BG_3)$coefficients[5,1], 
                        logLik(BG_3), with(summary(BG_3), 1-deviance/null.deviance)),digits=2)

TableA24[22,] <-round(c(nrow(BG3),coef(BG_3_cum)["pro_imms"],
                        coef(BG_3_cum)["pro_reds"],
                        coef(BG_3_cum)["nat_eus"], 
                        coef(BG_3_cum)["nat_regs"],
                        logLik(BG_3_cum), 1 - (logLik(BG_3_cum) / as.numeric(logLik(update(BG_3_cum, . ~ 1))))), digits = 2)
# SK
TableA24[23,] <-round(c(nrow(SK),summary(SK_3)$coefficients[2,1], 
                        summary(SK_3)$coefficients[3,1], 
                        summary(SK_3)$coefficients[4,1], 
                        summary(SK_3)$coefficients[5,1], 
                        logLik(SK_3), with(summary(SK_3), 1-deviance/null.deviance)),digits=2)

TableA24[24,] <-round(c(nrow(SK3),coef(SK_3_cum)["pro_imms"],
                        coef(SK_3_cum)["pro_reds"],
                        coef(SK_3_cum)["nat_eus"], 
                        coef(SK_3_cum)["nat_regs"],
                        logLik(SK_3_cum), 1 - (logLik(SK_3_cum) / as.numeric(logLik(update(SK_3_cum, . ~ 1))))), digits = 2)
# LV
TableA24[25,] <-round(c(nrow(LV),summary(LV_3)$coefficients[2,1], 
                        summary(LV_3)$coefficients[3,1], 
                        summary(LV_3)$coefficients[4,1], 
                        summary(LV_3)$coefficients[5,1], 
                        logLik(LV_3), with(summary(LV_3), 1-deviance/null.deviance)),digits=2)

TableA24[26,] <-round(c(nrow(LV3),coef(LV_3_cum)["pro_imms"],
                        coef(LV_3_cum)["pro_reds"],
                        coef(LV_3_cum)["nat_eus"], 
                        coef(LV_3_cum)["nat_regs"],
                        logLik(LV_3_cum), 1 - (logLik(LV_3_cum) / as.numeric(logLik(update(LV_3_cum, . ~ 1))))), digits = 2)
# EE
TableA24[27,] <-round(c(nrow(EE),summary(EE_3)$coefficients[2,1], 
                        summary(EE_3)$coefficients[3,1], 
                        summary(EE_3)$coefficients[4,1], 
                        summary(EE_3)$coefficients[5,1], 
                        logLik(EE_3), with(summary(EE_3), 1-deviance/null.deviance)),digits=2)

TableA24[28,] <-round(c(nrow(EE3),coef(EE_3_cum)["pro_imms"],
                        coef(EE_3_cum)["pro_reds"],
                        coef(EE_3_cum)["nat_eus"], 
                        coef(EE_3_cum)["nat_regs"],
                        logLik(EE_3_cum), 1 - (logLik(EE_3_cum) / as.numeric(logLik(update(EE_3_cum, . ~ 1))))), digits = 2)

TableA24


##++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# D.2 Comparing the Binary and the Adjacent Logistic Regression Models
##++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


#------------------------------------------------------------------------------------
# Table A26: Robustness check: Summary table. Binary and Adjacent Logistic Models
#------------------------------------------------------------------------------------

TableA26 <- matrix(NA, 28, ncol = 7)
colnames(TableA26) <- c("Obs.", "beta1", "beta2", "beta3", "beta4", "Log-Lik", "pse. R^2")
rownames(TableA26) <- c("Germany: binary", "Germany: adjacent", "France: binary", "France: adjacent", "Italy: binary", "Italy: adjacent", 
                        "Spain: binary", "Spain: adjacent", "Poland: binary", "Poland: adjacent", "Romania: binary", "Romania: adjacent", 
                        "Netherlands: binary", "Netherlands: adjacent", "Sweden: binary",  "Sweden: adjacent", "Hungary: binary", "Hungary: adjacent",
                        "Austria: binary", "Austria: adjacent", "Bulgaria: binary", "Bulgaria: adjacent", "Slovakia: binary", "Slovakia: adjacent",
                        "Latvia: binary", "Latvia: adjacent", "Estonia: binary", "Estonia: adjacent")

# DE
TableA26[1,] <-round(c(nrow(DE),summary(DE_3)$coefficients[2,1], 
                       summary(DE_3)$coefficients[3,1], 
                       summary(DE_3)$coefficients[4,1], 
                       summary(DE_3)$coefficients[5,1], 
                       logLik(DE_3), with(summary(DE_3), 1-deviance/null.deviance)),digits=2)

TableA26[2,] <-round(c(nrow(DE3),coef(DE_3_adj)["pro_imms"],
                       coef(DE_3_adj)["pro_reds"],
                       coef(DE_3_adj)["nat_eus"], 
                       coef(DE_3_adj)["nat_regs"],
                       logLik(DE_3_adj), 1 - (logLik(DE_3_adj) / as.numeric(logLik(update(DE_3_adj, . ~ 1))))), digits = 2)
# FR
TableA26[3,] <-round(c(nrow(FR),summary(FR_3)$coefficients[2,1], 
                       summary(FR_3)$coefficients[3,1], 
                       summary(FR_3)$coefficients[4,1], 
                       summary(FR_3)$coefficients[5,1], 
                       logLik(FR_3), with(summary(FR_3), 1-deviance/null.deviance)),digits=2)

TableA26[4,] <-round(c(nrow(FR3),coef(FR_3_adj)["pro_imms"],
                       coef(FR_3_adj)["pro_reds"],
                       coef(FR_3_adj)["nat_eus"], 
                       coef(FR_3_adj)["nat_regs"],
                       logLik(FR_3_adj), 1 - (logLik(FR_3_adj) / as.numeric(logLik(update(FR_3_adj, . ~ 1))))), digits = 2)
# IT
TableA26[5,] <-round(c(nrow(IT),summary(IT_3)$coefficients[2,1], 
                       summary(IT_3)$coefficients[3,1], 
                       summary(IT_3)$coefficients[4,1], 
                       summary(IT_3)$coefficients[5,1], 
                       logLik(IT_3), with(summary(IT_3), 1-deviance/null.deviance)),digits=2)

TableA26[6,] <-round(c(nrow(IT3),coef(IT_3_adj)["pro_imms"],
                       coef(IT_3_adj)["pro_reds"],
                       coef(IT_3_adj)["nat_eus"], 
                       coef(IT_3_adj)["nat_regs"],
                       logLik(IT_3_adj), 1 - (logLik(IT_3_adj) / as.numeric(logLik(update(IT_3_adj, . ~ 1))))), digits = 2)
# ES
TableA26[7,] <-round(c(nrow(ES),summary(ES_3)$coefficients[2,1], 
                       summary(ES_3)$coefficients[3,1], 
                       summary(ES_3)$coefficients[4,1], 
                       summary(ES_3)$coefficients[5,1], 
                       logLik(ES_3), with(summary(ES_3), 1-deviance/null.deviance)),digits=2)

TableA26[8,] <-round(c(nrow(ES3),coef(ES_3_adj)["pro_imms"],
                       coef(ES_3_adj)["pro_reds"],
                       coef(ES_3_adj)["nat_eus"], 
                       coef(ES_3_adj)["nat_regs"],
                       logLik(ES_3_adj), 1 - (logLik(ES_3_adj) / as.numeric(logLik(update(ES_3_adj, . ~ 1))))), digits = 2)
# PL
TableA26[9,] <-round(c(nrow(PL),summary(PL_3)$coefficients[2,1], 
                       summary(PL_3)$coefficients[3,1], 
                       summary(PL_3)$coefficients[4,1], 
                       summary(PL_3)$coefficients[5,1], 
                       logLik(PL_3), with(summary(PL_3), 1-deviance/null.deviance)),digits=2)

TableA26[10,] <-round(c(nrow(PL3),coef(PL_3_adj)["pro_imms"],
                        coef(PL_3_adj)["pro_reds"],
                        coef(PL_3_adj)["nat_eus"], 
                        coef(PL_3_adj)["nat_regs"],
                        logLik(PL_3_adj), 1 - (logLik(PL_3_adj) / as.numeric(logLik(update(PL_3_adj, . ~ 1))))), digits = 2)
# RO
TableA26[11,] <-round(c(nrow(RO),summary(RO_3)$coefficients[2,1], 
                        summary(RO_3)$coefficients[3,1], 
                        summary(RO_3)$coefficients[4,1], 
                        summary(RO_3)$coefficients[5,1], 
                        logLik(RO_3), with(summary(RO_3), 1-deviance/null.deviance)),digits=2)

TableA26[12,] <-round(c(nrow(RO3),coef(RO_3_adj)["pro_imms"],
                        coef(RO_3_adj)["pro_reds"],
                        coef(RO_3_adj)["nat_eus"], 
                        coef(RO_3_adj)["nat_regs"],
                        logLik(RO_3_adj), 1 - (logLik(RO_3_adj) / as.numeric(logLik(update(RO_3_adj, . ~ 1))))), digits = 2)
# NL
TableA26[13,] <-round(c(nrow(NL),summary(NL_3)$coefficients[2,1], 
                        summary(NL_3)$coefficients[3,1], 
                        summary(NL_3)$coefficients[4,1], 
                        summary(NL_3)$coefficients[5,1], 
                        logLik(NL_3), with(summary(NL_3), 1-deviance/null.deviance)),digits=2)

TableA26[14,] <-round(c(nrow(NL3),coef(NL_3_adj)["pro_imms"],
                        coef(NL_3_adj)["pro_reds"],
                        coef(NL_3_adj)["nat_eus"], 
                        coef(NL_3_adj)["nat_regs"],
                        logLik(NL_3_adj), 1 - (logLik(NL_3_adj) / as.numeric(logLik(update(NL_3_adj, . ~ 1))))), digits = 2)
# SE
TableA26[15,] <-round(c(nrow(SE),summary(SE_3)$coefficients[2,1],
                        summary(SE_3)$coefficients[3,1], 
                        summary(SE_3)$coefficients[4,1], 
                        summary(SE_3)$coefficients[5,1], 
                        logLik(SE_3), with(summary(SE_3), 1-deviance/null.deviance)),digits=2)

TableA26[16,] <-round(c(nrow(SE3),coef(SE_3_adj)["pro_imms"],
                        coef(SE_3_adj)["pro_reds"],
                        coef(SE_3_adj)["nat_eus"], 
                        coef(SE_3_adj)["nat_regs"],
                        logLik(SE_3_adj), 1 - (logLik(SE_3_adj) / as.numeric(logLik(update(SE_3_adj, . ~ 1))))), digits = 2)
# HU
TableA26[17,] <-round(c(nrow(HU),summary(HU_3)$coefficients[2,1], 
                        summary(HU_3)$coefficients[3,1], 
                        summary(HU_3)$coefficients[4,1], 
                        summary(HU_3)$coefficients[5,1], 
                        logLik(HU_3), with(summary(HU_3), 1-deviance/null.deviance)),digits=2)

TableA26[18,] <-round(c(nrow(HU3),coef(HU_3_adj)["pro_imms"],
                        coef(HU_3_adj)["pro_reds"],
                        coef(HU_3_adj)["nat_eus"], 
                        coef(HU_3_adj)["nat_regs"],
                        logLik(HU_3_adj), 1 - (logLik(HU_3_adj) / as.numeric(logLik(update(HU_3_adj, . ~ 1))))), digits = 2)
# AT
TableA26[19,] <-round(c(nrow(AT),summary(AT_3)$coefficients[2,1], 
                        summary(AT_3)$coefficients[3,1], 
                        summary(AT_3)$coefficients[4,1], 
                        summary(AT_3)$coefficients[5,1], 
                        logLik(AT_3), with(summary(AT_3), 1-deviance/null.deviance)),digits=2)

TableA26[20,] <-round(c(nrow(AT3),coef(AT_3_adj)["pro_imms"],
                        coef(AT_3_adj)["pro_reds"],
                        coef(AT_3_adj)["nat_eus"], 
                        coef(AT_3_adj)["nat_regs"],
                        logLik(AT_3_adj), 1 - (logLik(AT_3_adj) / as.numeric(logLik(update(AT_3_adj, . ~ 1))))), digits = 2)

# BG
TableA26[21,] <-round(c(nrow(BG),summary(BG_3)$coefficients[2,1], 
                        summary(BG_3)$coefficients[3,1], 
                        summary(BG_3)$coefficients[4,1], 
                        summary(BG_3)$coefficients[5,1], 
                        logLik(BG_3), with(summary(BG_3), 1-deviance/null.deviance)),digits=2)

TableA26[22,] <-round(c(nrow(BG3),coef(BG_3_adj)["pro_imms"],
                        coef(BG_3_adj)["pro_reds"],
                        coef(BG_3_adj)["nat_eus"], 
                        coef(BG_3_adj)["nat_regs"],
                        logLik(BG_3_adj), 1 - (logLik(BG_3_adj) / as.numeric(logLik(update(BG_3_adj, . ~ 1))))), digits = 2)

# SK
TableA26[23,] <-round(c(nrow(SK),summary(SK_3)$coefficients[2,1], 
                        summary(SK_3)$coefficients[3,1], 
                        summary(SK_3)$coefficients[4,1], 
                        summary(SK_3)$coefficients[5,1], 
                        logLik(SK_3), with(summary(SK_3), 1-deviance/null.deviance)),digits=2)

TableA26[24,] <-round(c(nrow(SK3),coef(SK_3_adj)["pro_imms"],
                        coef(SK_3_adj)["pro_reds"],
                        coef(SK_3_adj)["nat_eus"], 
                        coef(SK_3_adj)["nat_regs"],
                        logLik(SK_3_adj), 1 - (logLik(SK_3_adj) / as.numeric(logLik(update(SK_3_adj, . ~ 1))))), digits = 2)

# LV
TableA26[25,] <-round(c(nrow(LV),summary(LV_3)$coefficients[2,1], 
                        summary(LV_3)$coefficients[3,1], 
                        summary(LV_3)$coefficients[4,1], 
                        summary(LV_3)$coefficients[5,1], 
                        logLik(LV_3), with(summary(LV_3), 1-deviance/null.deviance)),digits=2)


TableA26[26,] <-round(c(nrow(LV3),coef(LV_3_adj)["pro_imms"],
                        coef(LV_3_adj)["pro_reds"],
                        coef(LV_3_adj)["nat_eus"], 
                        coef(LV_3_adj)["nat_regs"],
                        logLik(LV_3_adj), 1 - (logLik(LV_3_adj) / as.numeric(logLik(update(LV_3_adj, . ~ 1))))), digits = 2)

# EE
TableA26[27,] <-round(c(nrow(EE),summary(EE_3)$coefficients[2,1], 
                        summary(EE_3)$coefficients[3,1], 
                        summary(EE_3)$coefficients[4,1], 
                        summary(EE_3)$coefficients[5,1], 
                        logLik(EE_3), with(summary(EE_3), 1-deviance/null.deviance)),digits=2)

TableA26[28,] <-round(c(nrow(EE3),coef(EE_3_adj)["pro_imms"],
                        coef(EE_3_adj)["pro_reds"],
                        coef(EE_3_adj)["nat_eus"], 
                        coef(EE_3_adj)["nat_regs"],
                        logLik(EE_3_adj), 1 - (logLik(EE_3_adj) / as.numeric(logLik(update(EE_3_adj, . ~ 1))))), digits = 2)

TableA26

#-------------------------------------
# D.3 Redistribution versus Cohesion Policy
#-------------------------------------


#------------------------------------------------------------------------------------
# Table A28: Descriptive Statistics on Redistribution and Cohesion Policy
#------------------------------------------------------------------------------------

# Matrix
TableA28           <- matrix(NA, 14, ncol = 3)
colnames(TableA28) <- c("Redistribution","Cohesion Policy", "corr.")
rownames(TableA28) <- c("Germany", "France", "Italy", "Spain", "Poland", "Romania", "Netherlands", "Sweden", "Hungary", "Austria", "Bulgaria", "Slovakia", "Latvia", "Estonia")

# DE
TableA28[1,]  <-c(round(mean(DEycp$pro_red), digits=2),round(mean(DEycp$pro_cp), digits=2), round(cor(DEycp$pro_red,DEycp$pro_cp, method = "kendall"), digits=2))
# FR
TableA28[2,]  <-c(round(mean(FRycp$pro_red), digits=2),round(mean(FRycp$pro_cp), digits=2), round(cor(FRycp$pro_red,FRycp$pro_cp, method = "kendall"), digits=2))
# IT
TableA28[3,]  <-c(round(mean(ITycp$pro_red), digits=2),round(mean(ITycp$pro_cp), digits=2), round(cor(ITycp$pro_red,ITycp$pro_cp, method = "kendall"), digits=2))
# ES
TableA28[4,]  <-c(round(mean(ESycp$pro_red), digits=2),round(mean(ESycp$pro_cp), digits=2), round(cor(ESycp$pro_red,ESycp$pro_cp, method = "kendall"), digits=2))
# PL
TableA28[5,]  <-c(round(mean(PLycp$pro_red), digits=2),round(mean(PLycp$pro_cp), digits=2), round(cor(PLycp$pro_red,PLycp$pro_cp, method = "kendall"), digits=2))
# RO
TableA28[6,]  <-c(round(mean(ROycp$pro_red), digits=2),round(mean(ROycp$pro_cp), digits=2),round(cor(ROycp$pro_red,ROycp$pro_cp, method = "kendall"), digits=2))
# NL
TableA28[7,]  <-c(round(mean(NLycp$pro_red), digits=2),round(mean(NLycp$pro_cp), digits=2),round(cor(NLycp$pro_red,NLycp$pro_cp, method = "kendall"), digits=2))
# SE
TableA28[8,]  <-c(round(mean(SEycp$pro_red), digits=2),round(mean(SEycp$pro_cp), digits=2), round(cor(SEycp$pro_red,SEycp$pro_cp, method = "kendall"), digits=2))
# HU
TableA28[9,]  <-c(round(mean(HUycp$pro_red), digits=2),round(mean(HUycp$pro_cp), digits=2), round(cor(HUycp$pro_red,HUycp$pro_cp, method = "kendall"), digits=2))
# AT
TableA28[10,]  <-c(round(mean(ATycp$pro_red), digits=2),round(mean(ATycp$pro_cp), digits=2),round(cor(ATycp$pro_red,ATycp$pro_cp,method = "kendall"), digits=2))
# BG
TableA28[11,]  <-c(round(mean(BGycp$pro_red), digits=2),round(mean(BGycp$pro_cp), digits=2), round(cor(BGycp$pro_red,BGycp$pro_cp,method = "kendall"), digits=2))
# SK
TableA28[12,]  <-c(round(mean(SKycp$pro_red), digits=2),round(mean(SKycp$pro_cp), digits=2), round(cor(SKycp$pro_red,SKycp$pro_cp,method = "kendall"), digits=2))
# LV
TableA28[13,]  <-c(round(mean(LVycp$pro_red), digits=2),round(mean(LVycp$pro_cp), digits=2), round(cor(LVycp$pro_red,LVycp$pro_cp,method = "kendall"), digits=2))
# EE
TableA28[14,]  <-c(round(mean(EEycp$pro_red), digits=2),round(mean(EEycp$pro_cp), digits=2), round(cor(EEycp$pro_red,EEycp$pro_cp,method = "kendall"), digits=2))

TableA28
